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ABSTRACT 

We develop a theoretical model that explains the formation of hot coronae around 
strongly magnetized neutron stars — magnetars. The starquakes of a magnetar shear 
its external magnetic field, which becomes non-potential and is threaded by an electric 
current. Once twisted, the magnetosphere cannot untwist immediately because of its 
self-induction. The self-induction electric field lifts particles from the stellar surface, 
accelerates them, and initiates avalanches of pair creation in the magnetosphere. The 
created plasma corona maintains the electric current demanded by curl B and regulates 
the self-induction e.m.f. by screening. This corona persists in dynamic equilibrium: it 
is continually lost to the stellar surface on the light-crossing time ~ 10~^ s and replen- 
ished with new particles. In essence, the twisted magnetosphere acts as an accelerator 
that converts the toroidal field energy to particle kinetic energy. Using a direct nu- 
merical experiment, we show that the corona self-organizes quickly (on a millisecond 
timescale) into a quasi-steady state, with voltage e^e ~ 1 GeV along the magnetic 
lines. The voltage is maintained near threshold for ignition of pair production, and pair 
production occurs at a rate just enough to feed the electric current. The heating rate 
of the corona is ~ in agreement with the observed persistent, high-energy 

output of magnetars. We deduce that a static twist that is suddenly implanted into the 
magnetosphere will decay on a timescale of 1-10 yrs. The particles accelerated in the 
corona impact the solid crust, knock out protons, and regulate the column density of 
the hydrostatic atmosphere of the star. The transition layer between the atmosphere 
and the corona may be hot enough to create additional pairs. This layer is the likely 
source of the observed 100 keV emission from magnetars. The corona emits curvature 
radiation up to 10^^ Hz and can supply the observed IR-optical luminosity. We outline 
the implications of our results for the non-thermal radiation in other bands, and for the 
post-burst evolution of magnetars. 
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1. INTRODUCTION 

At least 10% of all neutron stars are born as magnetars, with magnetic fields B > 10^^ G. 
Their activity is powered by the decay of the ultrastrong field and lasts about 10^ years. They 
are observed at this active stage as either Soft Gamma Repeaters (SGRs) or Anomalous X-ray 
Pulsars (AXPs) (see Woods &: Thompson 2006 for a recent review). Twelve active magnetars 
are currently known, and all of them have similar parameters: spin period P ~ 10 s, persistent 
X-ray luminosity L ~ 10^^ — 10^^ erg/s, and surface temperature ^ 0.5 keV. The large-scale 
(dipole) component of the magnetic field i^dipoie > 10^^ G has been inferred from the measurement 
of rapid spindown in several AXPs, and more recently in two SGRs (Kouveliotou et al. 1998, 
1999). However, there is evidence for even stronger magnetic field inside the star, as was originally 
deduced from considerations of magnetic field transport (Thompson & Duncan 1996). In particular, 
the occurence of enormously energetic X-ray flares in SGRs points to internal magnetic fields 
approaching ~ 10^^ G (Hurley et al. 2005). Bursting activity has also been detected in two 
AXPs (Kaspi et al. 2003), which strongly supports the hypothesis that active rearrangement of the 
magnetic field takes place in AXPs as well as SGRs, although at differing rates. 

Besides the sporadic X-ray outbursts, a second, persistent, form of activity has been discovered 
by studying the emission spectra of magnetars. Until recently, the spectrum was known to have 
a thermal component with temperature k^T 0.5 keV, interpreted as blackbody emission from 
the star's surface. The soft X-ray spectrum also showed a tail at 2 — 10 keV with photon index 
r = 2 — 4. Intriguingly, the hardness of the 2 — 10 keV spectrum correlates with long-term 
activity as a transient burst source (Marsden &; White 2001). This deviation from pure surface 
emission already suggested that energy is partially released above the star's surface. Recently, 
observations by RXTE and INTEGRAL have revealed even more intriguing feature: magnetars are 
bright, persistent sources of 100 keV X-rays (Kuiper et al. 2004; Molkov et al. 2005; Mereghetti 
et al. 2005; den Hartog et al. 2006; Kuiper et al. 2006). This high-energy emission forms a 
separate component of the magnetar spectrum, distinct from the soft X-ray component that was 
known before. The new component becomes dominant above 10 keV, has a very hard spectrum, 
r ~ 1, and peaks above 100 keV where the spectrum is unknown. The luminosity in this band, 
L '-^ 10^^ erg s~^, even exceeds the thermal luminosity from the star's surface. These hard X-rays 
can be emitted only in the exterior of the star and demonstrate the presence of an active plasma 
corona. 

Our goal in this paper is to understand the origin and composition of the magnetar corona. The 
starting point of our analysis is the fact that an evolving magnetar experiences strong deformations 
of its crust, due to a shifting configuration of the internal magnetic field. These "starquakes" 
shear the magnetic field that is anchored in the crust, thereby injecting an electric current into the 
stellar magnetosphere. This process bears some resemblence to the formation of current-carrying 
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structures in the Solax corona, although there are important differences which shall be discussed 
below. We investigate how the starquakes and twisting of the magnetosphere can lead to the 
formation of a plasma corona around a magnetar. 

The persistence of magnctospheric twists can explain several observed properties of the AXPs 
and SGRs (Thompson, Lyutikov, & Kulkarni 2002; hereafter TLK). TLK investigated the obser- 
vational consequences using a static, force-free model, idealizing the magnetosphere as a globally 
twisted dipole. They showed that a twist affects the spindown rate of the neutron star: it causes the 
magnetosphere to flare out slightly from a pure dipole configuration, thereby increasing the braking 
torque acting on the star. The force-free configuration is independent of the plasma behavior in 
the magnetosphere. It does rely, however, on the presence of plasma that can conduct the required 
current = (c/47r)V x B. This requires a minimum density ric = Jb/gc. TLK showed that even 
this minimum density significantly modifies the stellar radiation by multiple resonant scattering, 
which can explain the tail of the surface emission at 2 — 10 keV. 

In this paper, we formulate the problem of plasma dynamics in the closed twisted magneto- 
sphere and find a simple solution to this problem, which allows one to understand the observed 
energy output of the magnetar corona and the evolution of magnetic twists. Plasma behavior 
around neutron stars has been studied for decades in the context of radio pulsars, and that prob- 
lem remains unsolved. The principle difference with radio pulsars is that their activity is caused by 
rotation, and dissipation occurs on open magnetic lines that connect the star to its light cylinder. 
Electron-positron pairs are then created on open field lines (e.g. Ruderman &; Sutherland 1975; 
Arons &; Scharlemann 1979; Harding & Muslimov 2002). By contrast, the formation of a corona 
around a magnetar does not depend on rotation of the neutron star. All observed magnetars are 
slow rotators (fJ = 27r/P ~ 1 Hz), and their rotational energy cannot feed the observed coronal 
emission. (The bolometric output of some magnetars is ~ 10^ times larger than the spindown 
luminosity, as deduced from the rotational period derivative!) The open field lines form a tiny 
fraction R^s^l/c ~ 3 x 10~^ of the magnetic fiux threading the neutron star. Practically the entire 
plasma corona is immersed in the closed magnetosphere and its heating must be caused by some 
form of dissipation on the closed field lines. In contrast with canonical radio pulsars, the magnetar 
corona admits a simple solution. 

Plasma near a neutron star is strongly constrained compared with many other forms of plasma 
encountered in laboratory and space, and therefore easier to model. The relative simplicity is due 

to the fact that the magnetic field is strong and fixed — its evolution timescale is much longer 
than the residence time of particles in the corona. Furthermore, particles are at the lowest Landau 
level and move strictly along the magnetic field lines, which makes the problem effectively one- 
dimensional. Finally, the magnetic lines under consideration are closed, with both ends anchored 
in the highly conducting stellar material. The basic questions that one would like to answer are 
as follows. How is the magnetosphere populated with plasma? The neutron star surface has a 
temperature k^T ^ 1 keV and the scale-height of its atmospheric layer (if it exists) is only a few 
cm. How is the plasma supplied above this layer and what type of particles populate the corona? 
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How is the corona heated, and what are the typical energies of the particles? If the corona conducts 

a current associated with a magnetic twist V x B 7^ 0, how rapid is the dissipation of this current, 
i.e. what is the lifetime of the twist? Does its decay imply the disappearance of the corona? 

An outline of the model proposed in this paper is as follows. The key agent in corona formation 
is an electric field E'y parallel to the magnetic field. E\\ is generated by the self-induction of the 
gradually decaying current and in essence measures the rate of the decay. It determines the heating 
rate of the corona via Joule dissipation. If E\\ = then the corona is not heated and, being in contact 
with the cool stellar surface, it will have to condense to a thin surface layer with k^T ^ 1 keV. The 
current-carrying particles cannot flow upward against gravity unless a force e£'|| drives it. On the 
other hand, when E^y exceeds a critical value, avalanches are triggered in the magnetosphere, and 
the created pairs screen the electric field. This leads to a "bottleneck" for the decay of a magnetic 
twist, which implies a slow decay. 

The maintenance of the corona and the slow decay of the magnetic twist are intimately related 
because both are governed by . In order to find i^y , one can use Gauss' law V • E = Airp, where 
p is the net charge density of the coronal plasma. This constraint implies that E and p must be 
found self-consistently. We formulate the problem of plasma dynamics in a twisted magnetic field 
and the self-consistent electric field, and investigate its basic properties. The problem turns out 
to be similar to the classical double-layer problem of plasma physics (Langmuir 1929) with a new 
essential ingredient: pair creation. 

We design a numerical experiment that simulates the behavior of the 1-D plasma in the mag- 
netosphere. The experiment shows how the plasma and electric field self-organize to maintain the 
time-average magnetospheric current j = demanded by the magnetic field, (47r/c) j = V x B. 
The electric current admits no steady state on short timescales and keeps fluctuating, producing 
avalanches. This state may be described as a self-organized criticality. Pair creation is found to 
provide a robust mechanism for limiting the voltage along the magnetic lines to e$e 1 GeV and 
regulate the observed luminosity of the corona. 

The paper is organized as follows. We start with a qualitative description of the mechanism of 
corona formation, and show that it can be thought of as an electric-circuit problem (§2). A careful 
formulation of the circuit problem and the set up of numerical experiment are described in § 3. 
In § 4 we find that the circuit without creation does not provide a self-consistent model. Pair 
creation is found to be inevitable and described in detail in § 5. In § 6 we discuss the transition layer 
between the corona and the relatively cold surface of the star. The maintenance of a dense plasma 
layer on the surface is addressed in § 7. Observational implications and further developments of 
the model are discussed in § 8. 
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2. MECHANISM OF CORONA FORMATION 

2.1. Ejection of Magnetic Helicity from a Neutron Star 

A tightly wound-up magnetic field is assumed to exist inside magnetars (e.g. Thompson &; 
Duncan 2001). The internal toroidal field can be much stronger than the external large-scale 
dipole component. ""^ The essence of magnetar activity is the transfer of magnetic helicity from the 
interior of the star to its exterior, where it dissipates. This involves rotational motions of the crust, 
which inevitably twist the external magnetosphere that is anchored to the stellar surface. The 
magnetosphere is probably twisted in a catastrophic way, when the internal magnetic field breaks 
the crust and rotates it. Such starquakes are associated with observed X-ray bursts (Thompson & 
Duncan 1995). The most interesting effect of a starquake for us here is the partial release of the 
winding of the internal magnetic field into the exterior, i.e., the injection of magnetic helicity into 
the magnetosphere. 

Since the magnetic field is energetically dominant in the magnetosphere, it must relax there 
to a force-free configuration with j x B = 0. The effective footpoints of the force-free field lines sit 
in the inner crust of the star. At this depth currents may flow across the magnetic lines because 
the deep crust can sustain a significant Ampere force j x B/c 7^ 0: this force is balanced by elastic 
forces of the deformed crust. The ability of the crust to sustain j -H B quickly decreases toward the 
surface, and it is nearly force-free where^ p ^ lO^'^iJ^g^ g/cm^ (Thompson &: Duncan 2001). The 
build-up of stresses in the crust from an initial magnetostatic equilibrium has been studied recently 
by Braithwaite & Spruit (2006). 

The ejection of helicity may be thought of as directing an interior current (which previously 
flowed across the field and applied a j x B/c force) into the force- free magnetosphere. The ejected 
current flows along the magnetic lines to the exterior of the star, reaches the top of the line and 
comes back to the star at the other footpoint. The currents emerging during the starquake may 
percolate through a network of fractures in the crust. Then large gradients in the field are initially 
present which can be quickly erased. A relatively smooth magnetospheric twist of scale ~ i^NS will 
have the longest lifetime. 

Globally twisted magnetic fields are observed to persist in the Solar corona over many Alfven 
crossing times. Twisted magnetic configurations have also been studied in laboratory experiments. 
In particular, the toroidal pinch configuration has been studied in detail and is well explained as 
a relaxed plasma state in which the imparted magnetic helicity has been nearly conserved (Taylor 
1986). The toroidal pinch relaxes to the minimum-energy state which has a uniform helicity density: 



^Existing calculations of the relaxation to magnetostatic equilibrium in a fluid star (e.g., a nasccut magnetar) 
assume that the initial toroidal and poloidal fluxes are comparable (Braithwaite & Spruit 2004). The initial toroidal 
flux can, in fact, be substantially stronger due to the presence of rapid differential rotation in the newly formed star. 

^Throughout this paper we use the notation X = X„ x 10", where quantity X is measured in c.g.s. units. 
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V X B = |LtB with jji = const. This state is achieved because the experiment has toroidal geometry 
and the current is isolated from an external conductor. The relaxation process is then free to 
redistribute the current between neighboring magnetic flux tubes via small-scale reconnection, 
which leads to a uniform /i. The geometry is different in the case of a magnetar (or the sun) — the 
field lines are anchored to the star and their configuration depends on the footpoint motion (for 
related calculations of the relaxed state of a driven plasma see, e.g., Tang Sz Boozer [2005]). The 
long lifetime 1 year) of a twisted configuration around magnetars is supported by observations 
of long-lived changes in the spindown rate and X-ray pulse profiles (Woods &; Thompson 2006; 
TLK) as well as persistent nonthermal X-ray emission. 

The shape of a twisted force-free magnetosphere may be calculated given the boundary condi- 
tions at the footpoints. An illustrative model of a globally twisted dipole field was considered by 
TLK. It has one parameter — twist angle Acp — and allows one to see the qualitative effects of 
increasing twist on the magnetic configuration. The current density in this configuration is 

j(r,^) = -^VxB~^sin2^A,^ (r > i?Ns), (1) 
47r 47rr 

where r, 9, (p are spherical coordinates with the polar axis chosen along the axis of the dipole. Even 
for strong twists, ~ 1, the poloidal components Br and Bg are found to be close to the normal 
dipole configuration, although the magnetosphere is somewhat inflated. The main difference is the 
appearance of a toroidal field B^. 

The global structure of a real twisted magnetosphere may be complex. The displacement ^ 
of the magnetosphcrc's footpoints during a starqiiakc may be described as a 2D incompressible 
motion with V • ^ = and V x ^ 7^ 0. Suppose that the yielding behavior of the crust produces a 
motion ^ 7^ that is localized within a plate S, and this plate is connected by magnetospheric field 
lines to another plate on the opposite side of the star. Then / V x ^ d^ = / ^ • dL = where L is 
the boundary of S. Hence V x ^ must change sign within S. For example, for an axially symmetric 
rotational motion, the ejection of current in the center of the plate is accompanied by initiation of 
an opposite current on its periphery that screens the ejected current from the unmoved footpoints 
with ^ = 0. 

In this paper, we study the plasma dynamics in a flux tube dS with a given twist current 

Jb = (c/47r)|V x B|. The magnetosphere may be thought of as a collection of such elementary flux 
tubes, whose currents js vary from footpoint to footpoint. The variation of Jb is expected to be 
gradual, so that the characteristic thickness of a flux tube with a given Jb is comparable to -Rns- 

The emerging currents are easily maintained during the X-ray outburst accompanying a star- 
quake. A dense, thermalized plasma is then present in the magnetosphere, which easily conducts 
the current. Plasma remains suspended for some time after the starquake because of the transient 
thermal afterglow with a high blackbody temperature, up to ~ 4 keV. The cyclotron resonance of 
the ions is at the energy heB/nipC = 6.3B15 keV and, during the afterglow, the radiative flux at 
the resonance is high enough to lift ions from the surface against gravity (Ibrahim et al. 2001). 



- 7- 



Eventually the afterglow extinguishes and the radiative flux becomes unable to support the 
plasma outside the star. The decreasing density then threats the capability of the magnetosphere to 
conduct the current of the magnetic twist. A minimal "corotation" charge density p^o = — fi- B/27rc 
is always maintained because of rotation of the star with frequency O (Goldreich & Julian 1969), 
but it is far smaller than the minimal density that is needed to supply the current j (eq. 1): 



Here Rmsix{f,0) = r/s\T?6 is the maximum radius of a dipolar field line that passes through 
coordinates (r, 0). The current, however, will be forced to continue to flow by its self-induction, 
and the magnetosphere finds a way to re-generate the plasma that can carry the current (see § 2.2). 
So, the twisted force-free configuration persists. 

The stored energy of non-potential (toroidal) magnetic field associated with the ejected current 
is subject to gradual dissipation. In our model, this dissipation feeds the observed activity of the 
corona. The stored energy is concentrated near the star and carried by closed magnetic lines with a 
maximum extension radius i?max ~ '^Rns- Thus, most of the twist energy will be released if these 
lines untwist, and we focus in this paper on the near magnetosphere -Rmax ~ "^Rns- 

2.2. Bottleneck for the Twist Decay and Plasma Supply to the Corona 

Consider a magnetic flux tube^ with cross section S ^ i?^g and length L ~ R^s which carries 
a current I = Sj. The stored magnetic energy of the current per unit length of the tube is 



The decay of this energy is associated with an electric field parallel to the magnetic lines Ey: this 
field can accelerate particles and convert the magnetic energy into plasma energy. Conservation of 
energy can be expressed as 



as follows from Maxwell's equations with E <^ B (e.g. Landau & Lifshitz 1975). The first term on 
the right-hand side is the Joule dissipation caused by Ey. The second term — the divergence of 
the Poynting flux — describes the redistribution of magnetic energy in space. 

The decay of the twist is related to the voltage between the footpoints of a magnetic line. 




(2) 




(4) 




(5) 



^In an axisymmetric magnetosphere, we could equally well focus on a flux surface of revolution which sits within 
a range AO of polar angle and has a cross section S = sin 0/S.O. 
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Fig. 1. — Schematic picture of a current-carrying closed magnetic flux tube, anchored to the star's 
surface. The current is initiated by a starquake that twists one (or both) footpoints of the tube. 
The current flows along the tube into the magnetosphere and returns to the star at the other 
footpoint, where it enters the atmosphere and the crust, following the magnetic field lines (j || B). 
Sufficiently deep in the crust, the current can flow across the fleld lines and change direction, so 
that it can flow back to the twisted footpoint through the magnetosphere (or directly through 
the star, depending on the geometry of the twist). A self-induction voltage is created along the 
tube between its footpoints, which accompanies the gradual decay of the current. The voltage is 
generated because the current has a tendency to become charge-starved above the atmospheric 
layer whose scale-height h = k^T/gnip is a few cm. 
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Here A and C are the anode and cathode footpoints and dl is the Hne element; the integral is 

taken along the magnetic line outside the star (Fig. 1). The product approximately represents 
the dissipation rate S in the tube. (The two quantities are not exactly equal in a time-dependent 
circuit.) The instantaneous dissipation rate is given by 

S 

I A 

where I(^) = 5j is the instantaneous current at position I in the tube. The current is fixed at the 
footpoints A and C, / = Iq; it may, however, fluctuate between the footpoints. 



rE-I(Odl, (6) 

J A 



" 3,. 



Note that the true scalar potential <I> is different from and not related to the dissipation 
rate. $ is the 0-component of the 4-potential = (<I>, A) which is related to the electric field by 

E=-V*-if. (7) 

If anode and cathode are connected by an ideal conductor then $(C) — ^{A) = while the voltage 
$e may not vanish. 

The voltage $e is entirely maintained by the self-induction that accompanies the gradual 
untwisting of the magnetic field. Thus, $e reflects the gradual decrease of the magnetic helicity in 
the tube, 

/ ABdV = - E-B5(/)d/. (8) 

J A J A 

(Here dV = S{l)dl is the volume element.) It also determines the effective resistivity of the tube: 
TZ = ^e/I- The self-induction voltage passes the released magnetic energy to charged particles, 
and a higher rate of untwisting implies a higher energy e$e gained per particle. A huge magnetic 
energy is stored in the twisted tube, and a quick untwisting would lead to extremely high Lorentz 
factors of the accelerated particles. 

There is, however, a bottleneck that prevents a fast decay of the twist: the tube responds to 
high voltages through the copious production of pairs. Runaway pair creation is ignited when 
electrons arc accelerated to a certain Lorentz factor ^± ~ 10^. (We describe the microphysics of 
pair creation in § 5.) The created plasma does not tolerate large — the plasma is immediately 
polarized, screening the electric held. This provides a negative feedback that limits the magnitude 
of $e and buffers the decay of the twist. 

A minimum $e is needed for the formation of a corona. Two mechanisms can supply plasma: 
(1) Ions and electrons can be lifted into the magnetosphere from the stellar surface. This requires 
a minimum voltage, 

e$e ~ gmRNS ~ 200 MeV, (9) 

corresponding to i^y that is strong enough to lift ions (of mass nii) from the anode footpoint and 
electrons from the cathode footpoint. (2) Pairs can be created in the magnetosphere if 

e$e ~ ^±mc'^ = 0.57± MeV, (10) 
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which can accelerate particles to the Lorentz factor 7-1- sufficient to ignite e creation. A third 
possible source — thermal production in a heated surface layer — will be discussed in § 6; the 
corresponding minimum voltage is currently unknown. 

If $e is too low and the plasma is not supplied, a flux tube is guaranteed to generate a stronger 
electric field. The current becomes slightly charge-starved: that is, the net density of free charges 
becomes smaller than |Vx B|/47re. The ultrastrong magnetic field, whose twist carries an enormous 
energy compared with the energy of the plasma, does not change and responds to the decrease of j 
by generating an electric field: a small reduction of the conduction current j induces a displacement 
current {l/4Tr)dEi/dt = (c/47r)V x B — j. The longitudinal electric field E^^ then quickly grows 
until it can pull particles away from the stellar surface and ignite pair creation. 

The limiting cases E'y (no decay of the twist) and E\\ — > 00 (fast decay) both imply a 
contradiction. The electric field and the plasma content of the corona must regulate each other 
toward a self-consistent state, and the gradual decay of the twist proceeds through a delicate 
balance: E^^ must be strong cnoTigh to supply plasma and maintain the current in the corona; 
however, if plasma is oversupplied E^<^ will be reduced by screening. A cyclic behavior is possible, 
in which plasma is periodically oversupplied and E^^ is screened. Our numerical experiment will 
show that such a cyclic behavior indeed takes place. 

3. ELECTRIC CIRCUIT 
3.1. Basic Properties 

Three facts facilitate our analysis of plasma dynamics in the electric circuit of magnetars: 

1. The ultrastrong magnetic field makes the particle dynamics 1-D. The rest-energy density 
of the coronal plasma is ~ mcj/e, where m is the mass of a plasma particle. It is much smaller 
than the magnetic energy density S^/Svr: their ratio is ~ mcP 

magnetic field lines are not perturbed significantly by the plasma inertia, and they can be thought 
of as fixed "rails" along which the particles move. The particle motion is confined to the lowest 
Landau level and is strictly parallel to the field. 

2. The particle motion is coUisionless in the magnetosphere. It is governed by two forces only: 
the component of gravity projected onto the magnetic field and a collective electric field E^^ which 
is determined by the charge density distribution and must be found self-consistently. 

3. The star possesses a dense and thin atmospheric layer. ^ Near the base of the atmosphere, 
the required current = (c/47r)V x B is easily conducted, with almost no electric field. Therefore 
the circuit has simple boundary conditions: E'y = and fixed current. The atmosphere is much 



''Such a layer initially exists on the surface of a young neutron star. Its maintenance is discussed in § 7. 
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thicker than the skin depth of the plasma and screens the magnetospheric electric field from the 
star. 

Our goal is to understand the plasma behavior above the screening layer, where the atmospheric 
density is exponentially reduced and an electric field E^^ must develop. The induced electric field 
E and conduction current j satisfy the Maxwell equation, 

Here j is parallel to the direction of the magnetic field. Projection of equation (11) onto the 
magnetic field gives 

^ = Mjb - j), JB = — |V X B|, (12) 

where we have taken into account the force-free condition V x B || B. If the conduction current j 
is smaller than js, then dE\\/dt > and an electric field appears that tends to increase the current. 
Alternatively, if j > Jb then an electric field of the opposite sign develops which tends to reduce 
the current. Thus, j is always regulated toward Jb- This is the standard self-induction effect. 

The timescale of this regulation, r, is very short. Consider a deviation \j — jsl Jb', then 
equation (12) gives 

E 

- 47rjB. (13) 

The current j = env (carried by particles with charge e and mass m) responds to the induced E'y 
by changing v, 

V eE» 



T m 



(14) 



Prom equations (13) and (14) we get 



~ = Up. (15) 



m 



Thus the timescale of the current response to deviations of j from Jb is simply the Langmuir plasma 
timescale cop^, which is very short. Taking the characteristic density of the corona, 

n^nc = — , (16) 
ec 



and Jb ~ (c/47r)(i3/r) (eq. 1), we estimate n ~ 10^^i?i5rg ^cm ^. The corresponding electron 



plasma frequency (m = me) is up ~ 10-*^^ Hz. The current relaxes to js on the timescale a;, ^ 



10 s, which is tiny compared with the light crossing timescale of the circuit. 

Local deviations from charge neutrality tend to be erased on the same plasma timescale. 
When the corotation charge density is neglected, the net charge imbalance of the corona must be 
extremely small, i.e. the corona is neutral. Indeed, the total number of pariclcs in the corona is 
Nc ~ ricR^s ~ 10^^ and a tiny imbalance of positive and negative charges 5Nc/Nc ~ 10~^^ would 
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create an electric field E ~ eSNc/R^g > GMrUp/eR^g that easily pulls out the missing charges 
from the surface. 

Adopting a characteristic velocity w ~ c of the coronal particles, the Debye length of the plasma 
Ad = v/ujp is essentially the plasma skin depth c/ujp. An important physical parameter of the 
corona is the ratio of the Debye length to the circuit size L ~ .Rns- 

In the case of the twisted dipolar magnetosphere (eq. [1]), this becomes 

C = 3 X 10-9 B'^^^^ {A(j) sin2 9) "^^^ . (18) 

The same parameter may be expressed using a dimensionless current Jb/J* ~ lO^^BisA^sin^ ^. 
Here j* corresponds to the plasma density n* = j*/ec such that Ad = L, 

m-r^ 7 -n:^ (19) 

Note that equation (15) for the plasma frequency is valid only for non-relativistic plasma. We will 
use the above definition of ujp, however, one should keep in mind that plasma oscillations in the 
corona may have a smaller frequency because of relativistic effects. 



3.2. Quasi-Neutral Steady State? 

Suppose for now that no are created in the corona, so that the circuit must be fed by 
charges lifted from the stellar surface. May one expect that a steady current j = Jb is maintained 
in the neutral corona by a steady electric force eE'y ~ grrii or e£'|| ~ ^mg? It turns out that such 
a balance between gravitational and electric forces is impossible in the circuit for the following 

reason: £"11 exerts opposite forces on positive and negative charges while the force of gravity has the 
same sign. This precludes a steady state where -Ey lifts particles into the magnetosphere without 
a significant deviation from neutrality appearing somewhere in the circuit. The point is illustrated 
by the following toy model. 

Suppose positive charges e (mass m_|_) fiow with velocity v+ > and negative charges — e (mass 
m_) flow with velocity v- < 0. The flow of each species is assumed to be cold in this illustration, 
with a vanishing velocity dispersion at each point, and the total current is j = j+ + j- where 
j± = ±en±v±. The continuity equation for +/— charges reads 

d(±en±/B) d(j±/B) ^ 

where I is the coordinate along the field line. Suppose the circuit is steady, d/dt = 0, and no strong 
deviation from neutrality occurs, i.e. n+ 2± n_ everywhere. Then one finds = —v+/v- = 
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const, that is the ratio of velocities must remain constant along the field line. Hence v+ and v_ 
must both vanish at the footpoints. This cannot be: for example, at the anode i^y must lift ions 
against gravity, and hence it is accelerating the electrons downward rather than stopping them^ 

So, if the cTirrcnt-carrying particles are supplied by the star, we expect either a significant 
charge separation at some places in the circuit, or a time-dependent behavior (or both). One must 
also allow for the creation of pairs if the light charges are accelerated to high enough energies, 
which will be discussed in § 5. Our approach to this problem is a direct numerical experiment that 
simulates the time-dependent behavior of plasma particles in the circuit. 



3.3. Numerical Simulation of a One-Dimensional Circuit 

We shall describe below a time-dependent simulation of a circuit operating along a thin mag- 
netic tube with given jB = (c/47r)|V x B| that is fixed in time. Strictly speaking, Jb is quasi-steady: 
it will be found to decay on a long resistive timescale yrs), much longer than the dynamic time 
of the current tdyn = L/c ^ 10~^ s. The purpose of the simulation is to understand how plasma is 
supplied above the star's surface layer, and to find the electric field that develops along the tube. 

The electric field may be decomposed as E = Ey + E^ where Ey is parallel to the magnetic 
field. It is E|| that governs the plasma motion and determines the release of magnetic energy in the 
circuit. In the fixed magnetic configuration, V x B does not vary with time and remains parallel 
to B (the force-free condition). Then equation (11) implies that dE/dt is parallel to B, i.e. no 
perpendicular field Ej^ is created by the induction effect, and V • Ej^ = may be assumed. Then 
Gauss' law reads 

47rp = V- E|| = (21) 

where I is length measured along the magnetic tube. The condition V • Ej^ = significantly 
simplifies the problem: it becomes strictly 1-D since Ex has no relation to charge density and falls 
out from the problem. Hereafter we assume E = Ey and drop the subscript ||. 

We note that the approximation of ID circuit, jsit) = const, excludes the excitation of 
transverse waves in the magnetosphere, which in reality can exist. These waves are described by 
the coupled fluctuations 6jB and E_l on scales much smaller than the circuit size L. (Wavelengths 
A ^ c/ujp can be excited by plasma oscillations.) These fluctuations may be expected to have a 
small effect on the circuit solution unless SJb becomes comparable to Jb- 

With jBit) = const, the particle motions on different magnetic lines are decoupled. Indeed, 
the particle dynamics is controlled by electric field E = E» which is related to the instantaneous 



®A static plasma configuration is possible with the electric and gravitational forces balancing each other at each 
point of the circuit. However, j = in this configuration and it is not relevant to our circuit problem. Besides, such 
a configuration would be unstable. 
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charge distribution by Gauss' law (21).^ Thus, E{1) on a given magnetic line is fully determined 
by charge density p{l) on the same line. The plasma and electric field evolve along the line as if 
the world were 1-D. 

The force-free condition x B together with V • = requires that j_b(/) oc B{1) along the 
magnetic line. We can scale out this variation simply dividing all local quantities (charge density, 
current density, and electric field) by B. This reduces the problem to an equivalent problem where 
— const. Furthermore, only forces along magnetic lines control the plasma dynamics, and the 
curvature of magnetic lines falls out from the problem. Therefore, we can set up the experiment 
so that plasma particles move along a straight line connecting anode and cathode (Fig. 2). We 
designated this line as the z-axis, so that I = z. 

This one-dimensional system may be simulated numerically. Consider N particles moving 
along the z-axis {N ~ 10^ in our simulations). Their positions Zi and velocities Vi evolve in time 
according to 



g{z) is the projection of the gravitational acceleration onto the magnetic line. In the simulations, 
we assume g = go{2z/L — 1) where L is the distance between the anode and cathode; then g = ±(70 
at the cathode/anode surfaces and g = in the middle z = L/2. The true gravitational acceleration 
depends on 7 and could be refined at large 7 by a factor ~ 2 using the exact geodesic equation. 
However, the crude Newtonian model is sufficient for our purposes. The electric field appearing in 
equation (22) is given by 



where E{0) = E{L) = are the boundary condition discussed above. Since the current j = Jb 
is fixed at both ends, the net charge between anode and cathode does not change with time and 
remains zero (in agreement with Given the instantaneous positions Zi of charges we immediately 
find E{z) from Gauss's law (eq. [23]). Thus, we have a well-defined 1-D problem of particle motion 
in a self-consistent electric field and the fixed gravitational field. To find a quasi-steady state of 
such a system we start with an initial state and let it relax by following the particle motion. 

The evolution of the system may be followed without calculating p or j: the equations contain 
only integrated charge Qiz) = pdz. In contrast with usual particle-in-ccll simulations, no grid or 
cells are needed. We find the exact E at any point z simply by counting the net charge of particles 
between z and anode. The only source of numerical error is caused by the finite timestep of particle 
dynamics, and we keep it small. At = 10~^ — 10~^ (L/c). 





(23) 



®The Maxwell equation (12) and Gauss' law (21) are equivalent when charge conservation dj/dl = —dp/dt is 
taken into account. 
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Fig. 2. — Set up of numerical experiment. Thin and dense plasma layers are maintained near the 
cathode and anode by injecting cold particles through the boundaries of the tube (the footpoints 
of a magnetic flux tube). The electric current is kept constant at the boundaries and the system is 
allowed to evolve in time until a quasi-steady state is reached. Without voltage between anode and 
cathode, the current cannot flow because gravity g traps the particles near the boundaries. The 
constant current at the boundaries, however, implies that voltage is immediately generated should 
the flow of charge stop in the tube. The experiment aims to find the induced voltage that keeps 
the current flowing. 
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The main dimensionless parameter of the problem is C = c/upL ~ 3 X 10-^ (see § 3.1). 
Although it is not possible to simulate circuits with such a low value of wc can experiment with 
a similar system with C, — 10^^. This allows one to understand the mechanism of the circuit and 
find plasma parameters that will be easy to scale to a real magnetosphere. A large number of 
particles N ^ 10^ may be followed with modern computers in a reasonable time, and C ^ 1 °iay 
be achieved by choosing an appropriate charge per particle. Two other basic requirements must 
be satisfied in the simulations: (1) the number of particles within the Debye length is much larger 
than one; and (2) the timestep is much smaller than uj'^^ and small enough to follow accurately 
the particle dynamics. Finally, we require the particle charge e to be sufficiently small, so that its 
electric field in our 1-D problem satisfies eE = 47re^ <C m(>g. Given a number of particles N , this 
last condition constrains C, from below. The minimum achieved in our numerical experiment is 
~ 5 X 10-=^. 

The ions and electrons in the simulated circuit have different masses rrii and mg. The ratio 
mi/rrie is smaller than it would be in a real magnetosphere. Below we will show models with 
mi/me = 10 and 30. Hereafter subscripts i and e refer to ions and electrons, respectively, and 
subscript p will refer to positrons. 

The neutral atmospheric layers near the boundaries of the computational box z = 0, L are 
created by injecting a flux of positive and negative charges with small temperatures Tg and Tj. 
Although most charges are bound to the star, some of them must be lifted by the electric field into 
the magnetosphere. The numerical simulation will show precisely how this happens. The boundary 
conditions E = Q and j = jb are imposed beneath the atmospheric layers. The condition j = jb 
implies keeping track of the flux of leaving charges, and injecting new particles with the correct 
imbalance between the positive and negative charges. Thus the required current is fed at a constant 
rate at the conducting boundaries. If net charge escaping through the cathode boundary during 
timestep dt is smaller than jsdt, we inject electrons to maintain the required jb, otherwise we inject 
the needed number of ions. Similarly, j = Jb is enforced at the anode boundary. This maintains 
the current, but does not determine the contributions to Jb from different types of particles. The 
plasma self-organizes during the experiment and decides itself what the contributions must be. 

The forced injection of net positive charge at anode and net negative charge at cathode would 
immediately create a long range electric fleld if the current were not flowing in the system. A 
deviation of j from js at any point of the circuit leads to appearance of E at that point. Thus, 
the boundary conditions j{0) = j{L) = js are equivalent to equation (12). 

The electric field does a net positive work on the current-carrying particles. In our experiment, 
the deposited kinetic energy sinks through the boundaries and never comes back. We envision 
dense radiativcly efficient layers just outside the simulation box and assume that the lost energy 
is emitted there. The surface layers maintained in the experiment have density only ~ 20 — 100 
times no = Jb/cc, which is just sufficient to screen the electric field and avoid interference of the 
boundaries in the circuit solution. A future full simulation may include radiative losses in situ in 



-17- 



the tube and allow one to track the fate of released energy and the spectrum of produced radiation. 

Our final choice regards the initial plasma state in the tube. We tried different initial states 
and got same results. In the runs shown below we take as initial conditions j{z) = js, a mildly 
relativistic bulk velocity Ve = Vi = 0.4c and Maxwellian momentum distribution with temperature 
equal to the bulk kinetic energy. 

The problem has four scales: tube length L, electron plasma scale Aoe = c/tope, ion plasma 
scale Adi = {mi/rrig)^^'^ \Df>, and scale-height of the atmospheric layer h = k{Te + Ti)/g(rni + nif,). 
For a real magnetosphere, Adc < Adi < h < L. We choose the parameters of the experiment so 
that these relations are satisfied. The inclusion of pair production in the simulation brings a new 
scale: mean free path I (see § 3.5 below). 

A numerical code implementing the described method was developed and first tested on simple 
plasma problems without gravity, g = 0. For example, we calculated the two-stream instability 
starting with two oppositely directed cold beams of positive and negative charges, and found the 
correct development of Langmuir oscillations on the plasma timescale. Various technical tests have 
been done, e.g. independence of the results of timcstcp and charge of individual particles e. The 
main parameter of the problem is and circTiits with equal and different N, e, dt gave the same 
result. Another test passed is the agreement of the linear-accelerator state obtained by the code 
(§ 5) with the solution of Carlqvist (1982). 

As we shall see below, creation is inevitable in the magnetospheric circuit. However, first 
we investigate what happens in the circuit with if pair creation "switched-off," when the current- 
carrying charges must be supplied by the atmospheric layers. 

4. CIRCUIT WITHOUT PAIR CREATION 

The circuit model has the following parameters: (1) dimensionless current js/j* = (see 
§ 3.1), (2) gravitational potential barrier (we will assume $g = c^/4 in numerical examples), 
(3) density of the hydrostatic boundary layer natm (which is regulated in the numerical experiment 
by the particle injection rate at the boundaries), and (4) injection temperatures of ions, Tj, and 
electrons, Tg. 

Note that Tj and Tg may be well above the surface temperature of the star k^T ^ 0.5 keV 
because the atmosphere is heated from above by the corona. The vertical temperature profile that 
is established below the corona is a separate problem which we do not attempt to solve in this 
paper. Tj and Tg serve as parameters of our model, and we explore first how the mechanism of the 
circuit depends on these parameters. 
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4.1. Thermally-Fed Circuit 

If the hydrostatic atmosphere is sufficiently hot, its Boltzmann tail at large heights may be 
able to conduct the current without the need to lift particles by an electric field. The plasma must 
remain quasi-neutral in the Boltzmann tail. One may think of the neutral atmosphere as a gas of 
composite particles m = mi + me with effective temperature T = Ti + T^. Its scale-height is given 

by 

h = = feB(?; + re) 

gom go{mi+me)' 

The plasma density is reduced exponentially on scale h. The reduced density remains high enough 
to conduct the current Jb all the way to the top of the gravitational barrier if 

natmexp(-0) >«c = ^. (25) 

One example of a circuit that satisfies condition (25) is shown in Figure 3. The electrons and 
ions are injected at the boundaries with fceTj = kBT^, = 0.5meC^ and mi = Wrrif,; this gives a 
scale-height h = 0.09L. After several light-crossing times L/c, the circuit develops the Boltzmann 
atmosphere that conducts the required current js = lO^j*. The electric field remains small and 
dynamically unimportant, eE < gouii. 

The thermally-fed regime requires too high a temperature for the pure ion-electron circuit. 
Indeed, mc^ = {rrii + me)c? ~ 1 GeV and then k^T > 1 GeV/ log(natm/?^c) would be required. 
Pair creation would take place at much lower temperatures (see § 6.3 below). Therefore, the 
thermally-fed pair-free circuit solution cannot apply to a real magnetosphere. The rest of this 
section focuses on circuits that do not satisfy condition (25). 



4.2. Linear Accelerator 



When the condition (25) is not satisfied, the circuit becomes plasma-starved and an electric 
field develops that lifts particles from the dense surface layer. We find that such a circuit quickly 
relaxes to a state where it acts as an ultra-relativistic linear accelerator. One example is shown 
Figure 4. The parameters of this particular model are js/j* = 10^, ksTg = 0.1^6?^ = 0.04meC^, 
and n-^tj^Jric — 20. We obtained similar solutions for various combinations of Tg, Ti, and natm/^^o 
that do not satisfy condition (25) of a thermally-fed regime. In all cases, the circuit reached the 
quasi-steady state at t ^ 5L/c. In this state, oscillations of electric field (on the plasma timescale 
LJp^) are confined to the thin atmospheric layers. A static accelerating electric field is created above 
the layers where the atmosphere density is exponentially reduced. The striking result is that on 
top of each layer, a large charge builds up, which is positive {+Q) at the anode and negative {—Q) 
at the cathode. The two charges create a long-range electric field E = AttQ. 
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Fig. 3. — Thermally-fed circuit without production. Upper panel: the plasma density n = ne+rii 
in units of ric = Jb/^c. Lower panel: deviation from neutrality p/en, and electric field in units of 
Eg defined by eEg = ruigo. In this simulation, = lOmg and C = c/ujpL ~ 0.03. The snapshot is 
taken at t = 20L/c, after the system has reached the quasi-steady state. 
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Fig. 4. — Circuit without e production. Upper panel shows the normahzed charge density and 
lower panel shows the electric field in units of Eg defined by eEg = rriigo. Solid and dotted curves 

correspond to two different moments of time, demonstrating that a quasi-steady state has been 
reached. In this simulation, mi = lOrue and (" = c/topL ~ 0.01. The final state is the relativistic 
double layer described analytically by Carlqvist (1982). 
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This configuration is a relativistic double layer (Carlqvist 1982). It is well described by Car- 

Iqvist's solution, which has no gravity in the circuit and assumes zero temperature at the bound- 
aries, so that particles are injected with zero velocity. According to this solution, the potential 
drop between anode and cathode, in the limit e$e ^ miC^, is given by 



where Z is the charge number of the ions (in our experiment Z = 1 was assumed). In the limit 
nii » Zm-e, this expression depends only on the ion and not the electron mass. 

The established voltage is much larger than is needed to overcome the gravitational barrier 

It does not even depend on as long as is large enough to prohibit the thermally-fed regime. 
Gravity causes the transition to the linear-accelerator state, but the state itself does not depend 
on ^g. The asymmetry of the double-layer solution (Fig. 4) is caused by the difference between the 
electron and ion masses. A similar but symmetric configuration is obtained for nii = rrie- 

The symmetric double layer with rrif. = nii is especially simple and can be qualitatively under- 
stood as follows. The currents carried by ions and electrons are equal by symmetry, ji = je = \j, 
and hence plasma is approximately neutral between the anode and cathode where Vi — c. A 
significant deviation from neutrality \ne — ni\ ^ ni + ng takes place in the acceleration regions Az 
near anode and cathode because here Vi differs from Ve- For instance, near the anode the ions are 
not accelerated yet to c, and their charge density pi = ji/vi is larger than that of the incoming 
electron beam, pg = je/ve — je/c, so the net charge density is positive. One can show that the 
anode charge Q{z) = jQ{ji/vi — je/ve)dz peaks at height z ~ Az where ions move with a mildly 
relativistic velocity Vi ~ c/2. Similarly, at cathode end, where the emitted electron flow is still 
slow, the net charge density is negative, and the charge —Q peaks at the distance Az from the 
cathode. The thickness Az and the characteristic electric field E in the mildly relativistic region 
are related by eEAz m^c? . Combining this with E = AttQ ~ AirencAz, one finds Az ~ Ad and 
Q ~ encAD- This is the only possible self-consistent static solution for a double layer (besides the 
trivial solution E = which is prohibited in our case by gravity). E = AttQ must be this large in 
a static double layer because a smaller E would accelerate the particles more slowly, which would 
imply larger non-neutral regions Az near the ends of the circuit and larger Q. Gauss' law would 
then imply a larger E, leading to inconsistency with the assumed smaller E. 

Similar estimates may be made for the asymmetric double layer with rrii > me- In the limit 
e$e S> rriiC?, one again finds ji ~ je from the Langmuir condition (see Carlqvist 1982), 



There is only one qualitative difference from the symmetric case. Now the mildly relativistic region 
near the anode {vi ~ c/2) has a thickness Aza ~ Adi, and near the cathode [ve ~ c/2) — 
Azc ~ Aoe- The spikes of charge near anode and cathode are now different: the cathode spike is 




(26) 




(27) 
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smaller by a factor of Adb/Adi = Z(nie/'miY^'^ . 
vanish to satisfy the boundary conditions E = 
between the two spikes, throughout the circuit, 
analytically by Carlqvist (1982). 



The net charge of a double layer, however, must 
0. The "missing" negative charge is distributed 
This behavior is observed in Fig. 4 and derived 



The static double-layer solution has been much studied previously, but its possible astrophysical 
applications remained unclear. Our numerical experiment shows that an ion-electron circuit in a 
gravitational field relaxes to the double layer of the macroscopic size L and huge voltage <I>e- The 
system does not find any state with a lower even though it is allowed to be time-dependent. The 
experiment shows that, regardless the initial state, the charges iQ near the anode/cathode build 
up with time in spite of the strong fluctuations that persist in the atmospheric layers (reversing 
the sign of E), and Q grows until the Carlqvist solution is reached. 

The effective boundaries of the double layer in our experiment are not exactly at 2; = 0, L 
because the atmosphere has a finite scale- height. Height z^, where the circuit becomes plasma- 
starved may be estimated from the condition riatm exp(— ~ Uc- 

If the linear accelerator were maintained in the twisted magnetosphere, the twist would be 
immediately killed off. The huge voltage implies a large untwisting rate and a fast dissipation of 
the toroidal magnetic energy (§ 2). The electron Lorentz factor developed in the linear accelerator 
is (taking the real rrii/me = 1836 and C ~ 3 x 10~^), 

^, = ^:^^^6xl0^ (28) 

However, new processes will become important before the particles could acquire such high energies: 
production of pairs will take place. Therefore, the linear-accelerator solution cannot describe 
a real magnetosphere. We conclude that pair creation is a key ingredient of the circuit that will 
regulate the voltage to a smaller value. 



5. PAIR CREATION IN THE CORONA 

In the magnetospheres of canonical radio pulsars with B ~ 10^^ G, pairs are created when 
seed electrons are accelerated to large Lorentz factors 7e ^ 10^. Such electrons emit curvature 
gamma rays that can convert to off the magnetic field. 

In stronger fields, another channel of creation appears. It is also two-step: an accelerated 
particle resonantly upscattcrs a thermal X-ray photon, which subsequently converts to a pair. This 
channel is dominant in radio pulsars with relatively strong magnetic fields (e.g. Hibschman & Arons 
2001). A similar channel of creation operates in the super-QED field of a magnetar. Because of 
the large magnetic field an accelerated electron can resonantly upscatter the ambient X-rays and 
produces pairs when its Lorentz factor is 7e ~ 10^, much smaller than in a normal pulsar. 
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5.1. Resonant Channel of Pair Creation 



Consider an electron moving along a magnetic line with Lorentz factor 7e 3> 1 and an ambient 
X-ray photon hcox propagating at an angle 6kB with respect to the electron velocity. In the electron 
frame, the photon is strongly aberrated and moves nearly parallel to the magnetic field with an 
energy huj'x = 7e(l — cos 9kB)fi^x- The resonant scattering may be viewed as a two-step process: 
the excitation of the electron into the first Landau level with energy 

1/2 



/ 2B 



rrieC 



^ i -Bqed 

followed rapidly by a de-excitation. Here Sqed = mic? /eh 4.4 x 10^^ G. 



(29) 



The resonant condition on photon energy hw'^ is obtained from energy and momentum con- 
servation, 

Ux + m,c^=^EB. ^ = ^pEE. (30) 



Here P is the velocity along B acquired by the excited electron (relative to its initial rest frame) 
and 7 = (1 — /3^)~^/^. The photon huj'x is assumed to propagate parallel to B, which is a very good 
approximation since we consider ultrarelativistic electrons and the ambient radiation is strongly 
beamed in the electron frame. Prom these equations one finds 7 = ^{EB/me(? + me(? /Eb) and 
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This resonant condition may also be written as a condition on 7e for given huJx and O^b, 



7r 



B/Bqed f hwx \ ^ 



1 - COS 6 kB KrUeC^ J 



lOkeV 



(31) 
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De-excitation is examined conveniently in the frame where the parallel momentum of the 
excited electron vanishes (Hcrold, Ruder, & Wunner 1982). The energy of dc-cxcitation photon E^ 
depends on its emission angle 6'em with respect to the magnetic field. This angle determines how 
the released energy Eb is shared between the emitted photon and the electron recoil. The relation 
E^{6em.) is found from energy and momentum conservation, 

E-y + {c^pI + mlc^f/'^ = Eb, E^ cos 6'em = "CPe- (33) 

Here Pe is recoil momentum of the electron, which is directed along B. Eliminating pe, one finds 

Eb 



E^{9em) 



sin^ Of. 



1 - 



1 4 

COS 6'em + ^§5- Sm^ ' 



^1 



1/2- 



In the super-critical field, B S> -Bqed; this simplifies to 

Eb 



Ej{9e 



1 -I- I COS ( 



(34) 



(35) 
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The de-excitation photon has the maximum energy = Eb when it is emitted at ^em = it/'2: (no 
recoil) and the minimum energy E^ = Eb/^ when ^em = (maximum recoil). 

Photons have two polarization states in the strongly magnetized background: 0-mode and 
E-mode. These eigenmodes are determined by the effect of vaccuum polarization in response to an 
electromagnetic wave. Both modes are linearly polarized: the 0-mode has its electric vector in the 
plane of the wave vector k and the background magnetic field B. The electric vector of the E-mode 
is perpendicular to this plane. De-excitation photons are emitted in either the O or E mode: the 
gyrational motion of the excited electron overlaps both linear polarization states. 

0-mode photons can convert to a pair with both particles in the lowest Landau state (Daugh- 
erty Sz Harding 1983). The threshold condition for this process is 

The threshold energy becomes 2me(? when transformed to the frame where the photon propagates 
perpendicular to B. Using equation (35), we find the interval of emission angles ^em for which the 
0-mode photon will immediately convert to a pair, 

/ 2\ 2 

l-|cOS^em| <8( -^j , 5>Bqed. (37) 

0-mode photons emitted below the threshold still have huge energy in the lab frame, ~ ^^-csEb, 
and will convert to e"^ after propagating a small distance Al ~ -RNs/7rcs- Indeed, even a photon 
emitted with ^em = will quickly build up a sufficient pitch angle for conversion because of the 
curvature of magnetic lines and/or gravitational bending of the photon trajectory. 

The threshold energy for an E-mode photon is -E^ in the frame where the photon propagates 
normally to the field. Therefore, an E-mode de-excitation photon cannot immediately convert to 
e^. It will, however, quickly split in two daughter photons, at least one of which will be in 0-mode 
(Adler et al. 1970) and will convert to a pair.^ 

Thus, resonant scattering in the ultrastrong field always leads to 6=*= creation, sometimes mul- 
tiple because E-mode may split repeatedly and give multiple 0-mode photons that convert to e^. 

The average multiplicity A4 is not much larger than unity, and M = I may be taken as a first 
approximation. The rate of pair creation by an accelerated electron is then determined by the cross 
section of resonant scattering and the number density of the target photons. The cross section of 
resonant absorption for a photon propagating parallel to B in the rest frame of the electron is given 
by (Daugherty & Ventura 1978), 



(38) 



^Splitting of the O-mode is kincmatically forbidden. All three components of the momentum are conserved in the 
splitting process. This allows splitting only for the mode with the smaller index of refraction, which is the El-mode 
in a super-QED magnetic field. 
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where re = e^/mgC^ and (TresO = 27r^e/5. This expression is valid for any B/Bqeb- 

Magnetic Compton scattering has also been calculated for arbitrary B/Bqed by Herold (1979) 
and Melrose &; Parle (1983). The calculations must be modified at the resonance to take into account 

the finite lifetime of the intermediate resonant state (e.g. Ventura 1979; Harding &: Daugherty 
1991). Using the results of Herold (1979) and Herold ct al. (1982), wc find the resonant Compton 
cross section that agrees with equation (38) within a factor of 2 (details will be given elsewhere). 
Below we use the cross section (38). Note also that at S » -Bqed, there are equal probabilities for 
the photon to end up in the E-mode or 0-mode in the final state, independent of ^em- 

The number density of target photons for resonant scattering is 

dLx /dlnujx ^ ^ ^ -,^.20 f ^Lx/dlnujx\ f fiujx\~^ _3 

nx ~ o « 2 X 10 — 5^ — — cm . (39) 

ATTR^^hwxc V 10^^ erg s"! J \keV J ^ ' 

The mean free path of an electron to scattering a photon is then 

^ in3 D f dLx/dlnuJx\~^ f f^x\ 

~ 10"" ^15 7:7:35 — -7— cm. (40) 

nxCTTesO V 10 ergs V \^eV J 

The upscattcrcd photons produce pairs at a small distance from the scattering site compared with 
Ires- So, the distance over which an accelerated electron creates a pair is approximately equal to 
the free path to scattering. It depends only on B and the spectrum of target photons. While being 
accelerated, an electron with initially low 7e may first experience a resonance with a high-energy 
photon. Therefore it is important to know the radiation spectrum at high energies, above the 
blackbody peak. It is possible that the transition region between the corona and the star emits 
100 keV photons (Thompson & Beloborodov 2005; see §§ 6 and 8.5 below), which may be targets 
for resonant scattering by electrons with relatively low 7e ~ 100. Photons with energy much above 
~ 100 keV are not good targets for resonant scattering in the strong-S region near the star because 
then Zres becomes comparable to i?NS- 



nxO'resO 



dlnB 



dlnl 



5.2. Inconsistency of the Linear Accelerator 

Now we can show that the linear accelerator discussed in § 4.2 becomes inconsistent in the 
presence of pair creation and cannot be established in the magnetosphere. Suppose it is established. 
Electrons in the accelerator will be in resonance with photons at the peak of the stellar blackbody 

spectrum, hwx ~ 1 keV, when they have Lorentz factors 7res ~ lO^-Bis. Using the Carlqvist 
solution for the electron Lorentz factor 7e(^), we find the characteristic distance from cathode 
where electrons reach 7res- This distance is given by 

Ores - TresAoe- (41) 

Here we have used Aoe ^ o-res ^ {frii/meY/'^L; the electric field of the double layer varies little 
across this region and equals E = 47r(j/c)ADe, where j is the current flowing through the circuit. 
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The resonance region is at the distance ares ~ 30 cm from cathode and its thickness is ~ Oj-es- The 
mean free path of an electron in this region is Zrcs ~ 10^ cm (sec eq. [40]), which is much larger 
than Ores- This means that only a small fraction / = Orcs/^rcs of the flowing electrons will upscatter 
a photon and the bulk of them pass through the resonance region without any interaction. 

In sum, / ^ 1 positrons are created per unit electron flowing through the circuit. The created 
positrons will flow to the cathode. This backflow carries a current fj and charge density p+ = fj/c. 
The charge density p+ creates a potential drop between the resonance region and cathode, 

3 

A$ = 27r^^. (42) 

'res ^ 

The linear-accelerator configuration becomes inconsistent if eA$ > ^j-^smeC? — then the created 
positrons will screen the assumed electric field. This condition may be written as 

lres<\lLCL. (43) 

We can evaluate it by substituting L = i?NS = 10 km and equation (40) for /res- Using the example 
of a twisted dipole magnetosphere with the corresponding expression (18) for one finds 

^ 7.. in32 R-V2^A^Al/2„^„fl ( 

din ojx 

The linear accelerator configuration is inconsistent if this condition is satisfied at some fuvx- It is 
easily satisfied at hujx ~ 1 keV (the observed luminosity in this range is ~ 10^^ erg/s). 



> 7 X 10^^ B-,'/' (A<^)V2 sin^ ' erg s-^ (44) 



Even if the linear-accelerator state were formally allowed as a self-consistent solution, the 
coronal circuit would not have to evolve to this state. Indeed, if the voltage in the circuit evolves 
gradually to higher and higher values, an discharge occurs well before the linear accelerator 
could be established. Wc now consider the circuit that operates on discharge, at a much lower 
voltage compared to the linear accelerator. 



5.3. Discharge via Pair Breakdown 

The development of pair breakdown may be illustrated with the following toy model. Suppose 
a uniform electric field E is fixed in the tube of size L. To start, we introduce one seed electron at 
the cathode with zero velocity. The electron accelerates toward the anode and its Lorentz factor 7e 
grows. The toy model assumes that the electron creates a pair when it reaches a certain threshold 
7res- At this point, the electron energy is shared between three particles, so the energy of the 
original electron is reduced. In this illustration wc use a random energy distribution between the 
three particles with the mean values of 1/2, 1/4, and 1/4 of ^yresfrieC^- This introduces stochasticity 
in the pair breakdown.^ 



In reality there are other sources of stochasticity, in particular, the random free-path to resonant scattering with a 
mean value Zres(7e)- It is taken into account in the numerical experiment described below (§ 5.3). For the illustrative 
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The key parameter is a-res/L where 

_ (7res - l)me(P 
^res — r-i 

eE 

is the length of electron (or positron) acceleration to the energy jj-esfneC^, starting from rest. The 
necessary condition for breakdown is ares/L < 1 — otherwise the seed electron reaches the anode 
before gaining enough energy for pair creation. 

The development of the pair breakdown is shown in a spacetime diagram in Figure 5. Each 

pair-creation event gives two new particles of opposite charge, which initially move in the same 
direction. One of them is accelerated by the electric field and lost at the boundary after a time 
< L/c, while the other is decelerated and can reverse direction before reaching the boundary. After 
the reversal, the particle is accelerated toward the opposite boundary and creates at least two new 
particles, one of which can reverse direction etc. This reversal of particles in the tube and repeated 
pair creation allows the plasma to be continually replenished. In the super-critical regime (left 
panel in Fig. 5) more than one reversing particle is created per passage time L/c and an avalanche 
develops exponentially on a timescale ~ fflrcs/c In the near-critical regime shown in the right panel 
of Fig. 5, just one reversing particle is created per passage time L/c. This critical state is unstable: 
sooner or later the avalanche will be extinguished. 

This toy model is very idealized because it assumes a fixed electric field E. In reality, the 
behavior of E is coupled to the charge density in the tube, which is determined by the particle 
dynamics and the boundary conditions j = Jb- E changes on a timescale much shorter than L/c 
and the behavior of the coupled plasma and electric field is complicated. 

The toy model, however, shows an essential property of the breakdown: it is a critical 
stochastic phenomenon. Above a critical voltage pair creation proceeds in a runway manner, and 
the current and the dissipation rate would run away if the voltage were fixed. Below the critical 

voltage pair creation does not ignite. The criticality parameter is L/a^^es = e<I'e/(7rcs — l)m(.c^. The 
tube with enforced current at the boundaries must self-organize to create pairs in the near-critical 
regime e$e ~ Tres^eC^ and maintain the current. 

If the current demanded by V x B is not maintained, the self-induction voltage will grow until 
it exceeds the critical value e$e ~ Tres'T^eC^ (corresponding to Ores ~ L) and a pair breakdown 
develops. Avalanches of then supply plasma that tends to screen the electric field by conducting 
the necessary current, and the voltage in the tube is reduced. The discharging tube is similar to 
other phenomena that show self-organized criticality, e.g. , a pile of sand on a table (Bak, Tang, & 
Weisenfeld 1987). If sand is steadily added, a quasi-steady state is established with a characteristic 
mean slope of the pile. The sand is lost (falls from the table) intermittently, through avalanches 
— a sort of "sand discharge." In oTir case, charges of the opposite signs are added steadily instead 



purposes of this section, it is sufBcient to introduce any random element in the pair creation process to make it 
stochastic. 
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of sand (fixed j at the boundaries), and voltage $e = EL plays the role of the mean slope of a 
pile. The behavior of the discharging system is expected to be time-dependent, with stochastic 
avalanches. 

A steady state is not expected for a pair-creating circuit for two reasons: (1) The created 
pairs may not maintain a static electric field. Like the Carlqvist double layer, a self-consistent 
static voltage must be huge, much above the critical value tolerated by pair discharges. No static 
solution exists at modest voltages because it implies a slow particle acceleration and a large charge 
imbalance in a broad region, leading to inconsistency (see the end of § 4.3). (2) The near-critical 
behavior is unstable: either a runaway is ignited or pair creation extinguishes soon after a seed 
particle is injected. 

Discharges are well known in tubes with low-temperature plasma, which is weakly ionized. 
In that case, the critical voltage accelerates electrons to the ionization threshold and avalanches 
of ionization develop. Charge carriers are then supplied by ionization rather than pair creation. 
Another difference is that the fast electrons are lost to the cylindrical wall of the tube, so the 
problem is not 1-D, and coherent waves of ionization can develop. 

5.4. Numerical Simulations 

5.4- 1- Setup 

We implement the pair-creation process in our numerical experiment in a simple way: electrons 
(or positrons) can create pairs if their Lorentz factors are in a specified range 7i < 7e < 72- The 
mean free-path for emission by an accelerated electron in this window, l^es, is a parameter of 
the experiment. It determines the probability of pair creation by the electron during a timestep dt: 
this probability is cdt/lj-es- Note that Adc ^ ^res ^ -Rns- 

In a real magnetosphere, the resonance window 71 < 7e < 72 is shaped by the spectrum of 
target photons and the local magnetic field B. In our simplified model, this window is constant 
along the magnetic line. It is possible to use a more realistic resonance condition that takes into 
account variations in B and F^^ along a magnetosphcric field line. However, our major goal is to 
see whether creation leads to a qualitatively different state of the circuit. We shall see that the 
new state weakly depends on the details of the resonance condition when B > Bqed and 7res 3> 1 
everywhere in the circuit. 

We consider here only circuits with boundary temperatures too low for plasma to be supplied 
thermally (§ 4.1). Then the circuit solution weakly depends on the chosen values of Tj and Tg. 
In this section, we show a sequence of models with ksTe = 0.04meC^ and fceTj = O.OAmiC^. The 
corresponding scale height of the surface layers is ^ = 0.04L. 

To illustrate the weak dependence on details of the pair creation process we consider two 
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extreme cases of very efficient and very inefficient e creation. Pair creation is less efficient if the 
resonant window (71,72) is narrow and the mean free-path l^es is large even in this window. We 
therefore focus on two models: Model A: an infinitely wide resonance window (7j-cs,oo) with free 
path ^rcs = and Model B: a narrow resonance window (7res) 27res) with a large mean free path 
^res = L/3. The inclusion of pair creation brings two new parameters into the experiment: 7res and 
Zres- We are especially interested in values of 7res that are comparable to or below mi /me, which 
corresponds to a real magnetosphere. 

54.2. Results 

The results of the numerical experiments are shown in Figures 6-14. Figure 6 displays the 
relation between current j and voltage <I>e in a circuit with rrii = lOmg and 71.cs = 10. The current 
is normalized to a minimum current j* for which the characteristic Debye length Adc = c/ujp equals 
the size of the system L (see eq. [19]). At small jb, the double layer configuration is established in 
the circuit; its voltage is not sufficient to ignite production. With increasing js, the voltage of 
the double layer increases as {jB/j*f'^'^ and reaches 7res'T^eC^ at ji ~ i*7res; then it stops growing. 
At jB > Ji the voltage saturates at the characteristic e$e ~ 7res^eC^ in both Models A and B. 
This asymptotic regime jb ^ Ji is of interest to us because magnetars are in this regime. All 
subsequent figures show the asymptotic behavior at jb ^ ji ■ 

Each point in Fig. 6 is obtained by running a time-dependent simulation that relaxes to a 
quasi-steady state. Figures 7 and 8 show the relaxation history in two simulations with js = 10"^ j^,. 
The voltage $e is zero in the initial state and grows on the dynamical timescale tdyn = L/c while 
the plasma and electric field self-organize to maintain the required current. The voltage grows until 
electrons in the circuit get accelerated to Lorentz factors 7e ~ 7res and pair production begins. One 
can see that after a few dynamical times the voltage stops growing and the circuit enters an 
oscillatory regime. During each oscillation, an increased voltage leads to a higher production 
rate iV+, then too many are produced which screen the voltage (discharge), drops, and 
begins to grow again. The oscillations persist until the end of the simulation (40-L/c), with a 
constant amplitude. 

We see that a quasi-steady state is reached after a few tdyn- This state is time-dependent on 
short timescales, but it has a well-defined steady voltage when averaged over a few tdyn- Models A 
and B have similar histories of $e and iV_|_ (cf. Figs. 7 and 8) except that Model A is noisier on 
short timescales and shows less coherent oscillations. (The quasi-periodic oscillations are, however, 
pronounced in the Fourier spectra of ^e[t\ and iV+[t].) It is not surprising that $e may vary on 
short timescales t <^ idyn- For example, if the current suddenly stopped, a characteristic voltage 
e$e = rUeC^ that is able to accelerate electrons would be created after time t = nieC^ /^T^eLjB = 
{jB/j.)-\L/c)^L/c. 
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In all circuits of interest (jb » ji) we found the time-average voltage, 

e$e ~ Tres^^eC^- (45) 

This relation applies even to circuits with 7resWeC^ S> w-jC^, where lifting of the ions is energetically 
preferable to pair creation. We conclude that the voltage along the magnetic tube is self-regulated to 
a value just enough to maintain pair production and feed the current with pairs. The robustness 
of this result is illustrated in Figure 9, which shows circuits with mi/m.(, = 10,30 and varioTis 7ics- 
In all cases, e$e ~ IvcsT^eC^ is established, and the pair creation rate 2A^+ ~ Jb/c is maintained. 
Thus, the near-critical state is self-organized as expected. Avalanches of creation are observed 
to happen anywhere in the circuit between the dense boundary layers. 

The current is carried largely by everywhere in the tube. The ion fraction in the current 
depends on the ratio e^^e/fnic^ ~ jres'me/nT'i (Fig- 10). If this ratio is small, pairs arc easily 
produced with a small electric field, too small to lift the ions, and the ion current is suppressed. 
In the opposite case, ions carry about 1/2 of the current. For typical parameters expected in the 
magnetosphere, 7res ^ mi/me, ions carry ^ 10% of the current. 

The simulations show that the released energy sinks through the anode and cathode boundaries 
at approximately equal rates, and that their sum equals js^e- For magnetars, this implies that 
both footpoints of a twisted magnetic tube will radiate energy bolometrically at comparable rates. 



5.4-3. Sample Model 

We now describe in more detail the plasma behavior in one sample experiment. It assumes 
Model A for creation and 7i.cs = mi/rrie = 10.^ The relaxation history of this experiment is 
shown in Figure 7, and a snapshot of the circuit is shown in Figure 11. 

The charge density p = e{ni + rip ~ n^,) and electric field arc found to fluctuate significantly in 
time, while the densities n^, n^, and are approximately steady. Deviations from charge neutrality 
are small: p / e{ni + Up + Ui) <C 1 and E'^/Snllpi <^ 1 where Upi is the plasma kinetic energy density. 
The plasma turbulence may then be described as a superposition of Langmuir waves. We do not, 
however, attempt a development of analytical theory in this paper. 

The anode current je + ji + jp is dominated by electrons. The anode value of je even exceeds 
js by 17%. Together with ions {ji ~ O.lSjs), electrons balance the negative current of positrons 

jp ~ —O.Sjs so that the net current je + ji + jp equals js- Positrons contribute negative jp at the 
anode because they move in the "wrong" direction (they are created by accelerated electrons, and 
some of them leave the box against the mean electric field). The cathode current is dominated by 
positrons (64%); electrons contribTitc 23%, and ions 13%. 



®A similar solution is obtained for 7res = mi /me = 30. We expect the same circuit solution to apply to 7res = 
mi/me = 1836 or any other 7res = mi/me ^ 1. 
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Bulk speeds of the three plasma components, Vg, Vi, Vp, and the dispersions of their dimen- 

sionless momenta p/mc = 7/? are shown in Figures 12 and 13. The dispersions are significant, i.e. 
the corresponding distribution functions arc broad. The timc-avcragcd distribution of pair creation 
over z is shown in Figure 14. It has two peaks that correspond to pair production by accelerated 
e~ and e"^. The net rate of pair creation by e" is about 2 times higher than that by e"^ in this 
model. 

A complicated kinetic state self-organizes in the phase space of electrons, ions, and positrons. 
The state keeps oscillating. The oscillations might be related to the fact that positrons move in 
a wrong direction at creation (negative current) and have to reverse direction to maintain the net 
j = Jb without strong deviations from neutrality. This requires an extra voltage which leads to 
"overshooting" in the pair creation rate followed by a temporary screening of the electric field. 



6. TRANSITION LAYER 

A transition layer must be established between the corona and the relatively cold optically 
thick atmosphere of the star. It may play a crucial role for two reasons. (1) The layer is able 
to radiate away the energy received from the corona because it is dense enough for efficient 
bremsstrahlung emission. The corona itself may not efficiently radiate the dissipated power be- 
cause its bremsstrahlung output is negligible, and Compton cooling is inefficient as upscattered 
photons cannot escape — they convert to pairs in the ultra-strong magnetic field. (2) The layer 
may be hot enough to maintain thermal pair creation at an interesting rate and feed the charge 
flow in the coronal circuit (Thompson & Beloborodov 2005). A detailed analysis of the transition 
layer is left for future work. Below we describe some key mechanisms by which energy is deposited 
and redistributed within the layer, and by which pairs may be created. 



6.1. Heating by Beam Instability 

The particles accelerated in the corona enter the hydrostatic atmosphere with high Lorentz 
factors 7f, ~ e^e/iT^eC^ ■ This relativistic beam drives an electrostatic instability with the growth 
length (e.g. Godfrey, Shanahan, & Thode 1975), 

l^r^jJ !h] ^ = 1^ ( . (46) 

\nbj ujpe (47rn6re)V2 \nej 

Here LUpg = (Aimee^ /meY^"^ is the electron plasma frequency in the atmosphere with density ng, 
and Te = e^/meC^ — 2.8 x 10~^^ cm. The beam exponentially amplifies a seed Langmuir wave whose 
phase speed is in resonance with the beam. The exchange of energy with the atmospheric plasma 
occurs on a scale comparable to l^. As the beam propagates a distance ~ into the layer, the 
amplified Langmuir waves feed back on the beam and flatten its distribution function to a plateau. 
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If the initial distribution function of the beam is ideahzed by (5-function, fbile) = const 5{je — Tb); 
its relaxation to the plateau 1 < 7e < 75 implies the deposition of ^ of the beam energy into 
plasma turbulence. Since the relaxation occurs in a strong magnetic field, the particle momentum 
distribution remains one-dimensional. This is different from beam relaxation in a normal plasma 
— in that case a large dispersion of perpendicular momenta develops before relaxing to the plateau 
in energy. 

The coronal beam has a density n;, ~ ric = js/ec ~ 10-^^ — 10^^ cm"^ and Lorentz factor 75 ~ 
7res ~ 10^ (see § 5; this would correspond to a total energy deposition rate of 10^^ 76,3?^i6 erg 
if integrated over the entire surface of the star). Substituting these values to equation (46) one 
finds lb ^ ^ cm. The beam relaxation occurs at a much larger depth than the skin depth of the 
atmospheric plasma, Ih ^ c/ujp. Therefore, the beam instability may heat the atmosphere only 
well below the screening layer, i.e. well below the region where the induced electric field develops 
and particles are lifted to the corona. 

In our numerical experiments ~ 10(c/a;p) ~ O.IL was larger than the scale height of the 
atmosphere, h, so heating by beam instability could not possibly operate in the computatonal box. 
The possible atmospheric heating at larger depths was regulated in the experiment by hand, by 
setting the boundary temperatures Tg and Tj. Further investigation of beam heating is deferred to 
a future work. 

Comparing the expected li, and /i, one has 



Note that l^/h is reduced if the plasma temperature grows. It is possible that the atmosphere 
develops a hot layer with 1^ -C h, which absorbs a significant part of the beam energy. The excited 
Langmuir turbulence heats mostly electrons and may grow to k^Te > 100 keV. Then hard X-ray 
emission and thermal pair creation may occur as discussed below. 



Dissipation without losses of energy would imply an unlimited growth of temperature. If 
the twisted magnetosphere were thermally insulated, it would quickly heat up and remain filled 
with a hot plasma that would steadily conduct the current with no further dissipation. Losses 
of the dissipated energy are therefore crucial in the regulation of the state of the circuit. In our 
numerical experiment we simply allowed the released energy to sink through the anode and cathode 
boundaries, assuming that it is radiated away outside our boundaries. Here we discuss a possible 
mechanism of the energy loss. 

We expect the energy of the accelerated coronal particles to be deposited in an optically thin 
atmospheric layer by collisionless processes (such as beam instability discussed in § 6.1). This layer 




(47) 



6.2. 



Thermal Conduction and Radiative Losses 
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is still unable to cool radiatively, and the deposited heat will be conducted to deeper and denser 
layers that can radiate the heat away by bremsstrahlung (Thompson & Bcloborodov 2005). The 
radiation flux emitted at that depth is ^ heg where /i is a characteristic scale-height and eg 
is the free-free emissivity. Plasma with temperature Tg and density rig has the bremsstrahlung 
emissivity 

es^aey^ila^m.c^ 9^^'^, (48) 

where a = /he = 1/137 is the fine structure constant. The characteristic scale-height on which 
density and temperature change in the layer is the hydrostatic scale-height, h = kBTg/mig (we 
assume here Tg < Tj since the collisionless dissipation heats mostly electrons). Then the emitted 
radiation flux is 

2 

F^ = heg^a — d^/'^—nla^meC^. (49) 

rrii g 

It carries away a significant fraction of the energy flux received from the corona when F^ ^ F = 
ficlbfneC^, where 7b = e^ej^e^? — 7res- This gives one relation between temperature Qe and 
Thomson optical depth tt = ha^ne at which the heat flux is converted to radiation, 

0-^/2,2 ^-JLl^Hh^^ (50) 

where 

Tc = -RNSf^xnc ~ 0.01 - 0.1. (51) 

Equation (50) is not sufficient to determine Be and tt of the emitting layer because the mech- 
anism of heat transport to this layer has not been specified yet. A second relation is needed, which 
is the equation of thermal conductivity (Thompson & Bcloborodov 2005). It relates the heat flux 
Fc ^ F F^ to the temperature gradient: Fc = ndTc/dz where k = UeVcikB, Vc = (/ceTe/me)^/^, 
and i is the mean free path of electrons. In the absence of plasma turbulence, the mean free path 
is determined by Coulomb collisions, 

^Coui = (52) 

rieCTT 

If strong plasma turbulence develops with energy density above the thermal level, Utuvh = {^F^^)'^/8Tr ^ 
^BTe/A^jg, the fluctuating electric field suppresses the mean free path. It is then given by 

^^ADef^l^V'. (53) 

\nekBTe ) 

(We here assumed that the turbulence develops at frequencies uo ~ wpe.) We will keep ^ as a 
parameter in our estimates. It is bounded from above and below by Aoe < i < ^Coui- 

Estimating dTg/dz ~ Tg/h, one finds 

Fc = HeVei—r— =neVei guii, (54) 
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- - ij^) . (55) 

Combining equations (50) and (55), one finds the temperature of the radiating layer, 

'^''['^^^^^^.-r) ' ^BTe. 200 (-) (^) [j-^ keV. (56) 

If ^ < ^Coui (suppressing heat conduction) the emission temperature will grow above 200 keV. The 
growth, however, will be quickly stopped at ~ 1 because the hot layer starts to generate pairs 
which can short out the circuit (see § 6.3). 



6.3. Heated Atmospheric Layer as a Possible Source of e 

6.3.1. Channels of Pair Creation 

The dominant channels of thermal pair creation in a hot atmosphere involve an electron-ion 
collision. Three such channels are summarized below. 

1. An pair can be created directly in a collision between an ion and an energetic electron (or 
positron), e" + Z ^ e~ + Z + e~ + e'^ . This reaction is possible when the kinetic energy of the 
incident electron exceeds a threshold £?thr = 2meC^ (both e~ and e"^ can be created in the lowest 
Landau state). When the temperature k^Te ^ 2meC^, the reaction occurs in the high-energy tail 
of the electron distribution. 

2. A bremsstrahlung photon can be emitted in an electron-ion collision, almost always in the O- 
mode. The photon converts to an e"^ pair if satisfies the threshold condition > 2me(? / sin^em 
where ^em is the emission angle (cf. § 5.1). 

This second channel of pair creation is similar to the first one, where a virtual photon converts 
directly to a pair during the clcctron-ion collision. The emission of a real bremsstrahlung photon 
makes the process two-step and increases its threshold from 2meC^ to 2meC^/ sin6'em- The exact 
cross sections of these processes in the ultrastrong magnetic field have not been calculated, to our 
knowledge. They may be roughly estimated using the nonmagnetic formulae (e.g. Berestetskii, 
Lifshitz, &; Pitaevskii 1982): (Jbremss ~ Q.SZ^aaT and (Jdirect ~ QiCbremss- 

Note that relativistic bremsstrahlung emission is preferentially beamed along B, in the direction 
of the incident electron, so ^em is small. Then the threshold is substantially higher than 2me(? and 
the bremsstrahlung photon will be unable to convert to in situ. If the photon propagates away 
from the emission site where B has a different direction, the threshold may be reduced because 
the pitch angle of the photon increases. Gravitational bending of the photon trajectory also helps 
to increase the pitch angle. Thus, pair creation by bremsstrahlung photons can be non-local and 
take place above the hydrostatic atmosphere. 
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3. An electron can be excited to the first Landau state by a Coulomb collision with an ion. The 
excited electron rapidly de-excites and emits an energetic photon that may convert to a pair (or 
split, see § 5.1). This channel requires the initial electron to have a minimum energy E'thr = Eb ~ 
6.75^5^ mgC^ (eq. 29). The side-scattering into the first Landau state involves a perpendicular 
momentum transfer p± = {2B / Bq^DY/'^nieC, and the corresponding Coulomb cross section may 
be estimated as 

-cou.(0 - 1) - = — (^^J ar. (57) 



6.3.2. Feeding the Corona with Atmospheric Pairs 

The rate of pair creation in the atmosphere may be roughly estimated as the rate of electron-ion 
collisions above the threshold -Ethr, 

n+ ~ Ui /thr Ue a c, (58) 

where /thr ~ exp(— E'thr/^B^e) is the fraction of thermal 1-D electrons that have energy > -Ethr; 
(7 and £^thr are the cross section and the threshold energy of the relevant channel. Atmospheric 
pair production is able to feed the coronal current jb if hh^ > Jb/s ot, equivalently, 

/thr TTUi j > "c, (59) 

where tt = a^Ueh. This condition may be rewritten using Tc = -Rns^t^^c and h = 6e{me/mi){(? / g), 

rrip 



^e-V| > ^ f-^^) ^ ^ - 2 X 10-3 ^ ^. (60) 

rrip \RmgJ O- /thr /thr o- 



If this condition is met at some depth tt, the atmospheric pairs can dominate the charge supply 
to the coronal circuit. At tt » 1 the parameter OgT^ ^ 1 is limited by Compton cooling of the 
plasma, which implies a low Tg and negligible The condition (60) may be satisfied only at 

Tt ^ 1 and requires at least mildly relativistic Tg. The created pairs have kinetic energies ~ fce^e) 
so they may elevate, overflow the gravitational barrier, and maintain the current. 



6.3.3. Numerical Models 

In order to sec how atmospheric pair production could affect the circuit, we have done a few 
numerical simulations using the following toy model. We assume that the atmosphere responds to 
the heat flux from the corona by producing pairs. Therefore, in addition to the ion-electron 
atmosphere, we now maintain a thermal atmosphere. Its density at each boundary is regulated 
in the experiment so that it is proportional to the flux of lost energy through the boundary. 

If the atmospheric pairs were cold, the circuit would be in the same regime as in their absence. 
The fact that positrons are now available at the boundaries to carry the current instead of ions does 
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not lead to any significant reduction of the voltage: it still evolves to the state with e$e — 7res"^eC^5 
in which energetic are supplied via the resonant-scattering channel. We have checked this by 
running the models with cold atmospheres. The key property of atmospheric pairs is, however, 
that they are created hot enough to feed the current thermally (cf. § 4.1). Their expected kinetic 
energy is comparable to the height of the gravitational barrier (§ 6.3.2) and they can overflow the 
barrier. Therefore, the required current Jb can be conducted at a small voltage e$e < Tres^JeC^- 

The thermally-driven regime was indeed established in our experiment with a realistic /cb7± = 
O.bmeC^- The established voltage then depends on the efficiency of pair supply in response to the 
energy flux from the corona. We deflne the efficiency ^ as the fraction of the energy flux through 
the boundary that is returned to the corona in the form of kinetic and rest-energy. We find that 
the voltage $e across the circuit adjusts so that the flux of released energy, Jb^b, can generate the 
minimum flux of atmospheric pairs to feed the coronal current, Jb/^- Then voltage is related to 
^ roughly as e$e ~ ^"^mgC^. For example, a maximum efficiency ^ ~ 1 would reduce the voltage 
to e$e ^ rUeC^- We have done the corresponding experiment that confirmed this estimate. The 
voltage was found to fliictTiate significantly, but it did not get high enough to initiate the resonant 
channel of production (7res = rrii/me = 10 in the simulation). The ion current was small, 

ji/jB ^ 10-2. 

When ^ is reduced, a higher voltage is required in the circuit according to the relation e$e ~ 
^-^mgC^. This estimate is applicable as long as ^ ^ Trcs' that e<I>e < Ixcsi^eC? ■ A further 
reduction of ^ does not lead to an increase of voltage: then the circuit is mainly fed by pairs 
created via the resonant channel, and e$e ~ Tres^ieC^ as described in § 5. The thermal pair 
creation in the atmosphere then has a negligible effect on the circuit. 

We have also studied one modification of this experiment by allowing only the anode to supply 
pairs with efficiency ^ = 1/2. In this case, the atmosphere created by the anode filled the whole 
tube, and conducted the current. Ions lifted from the anode surface contributed a small fraction to 
the current ji/jB — 3 x 10-^. The established voltage was ~ 5 times higher than in the model with 
pair supply at both ends. However, the circuit was able to operate without igniting the resonant 
pair creation (which we checked by setting 7res = oo)- This shows that the thermally-driven regime 
can operate even when pair creation takes place at one footpoint of the magnetic tube. 

In summary, the coronal current may be fed by atmospheric instead of the pair discharge 
via the resonant channel if ^ ^ 7res- T^^^ numerical coefficient in this condition increases if one of 
the surfaces (anode or cathode) dominates the pair supply. 

6.4. Anomalous Resistivity and Hydrostatic Equilibrium 

In addition to the heating by the relativistic beam that comes from the corona, the atmosphere 
may experience anomalous heating in response to the strong electric current. This type of heating is 
known to occur in astrophysical and laboratory plasmas when the drift speed of the current-carrying 
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electrons, Vd, exceeds the sound speed Cg = {k^Tf./miY/'^ (e.g. Bernstein & Kulsrud 1961). If the 
electron current can be described as a shift of a MaxwelHan distribution with respect to the ions 
in velocity space, A^e = Vd, a positive gradient dfc/dv^, > appears in the electron distribution 
function viewed from the rest-frame of the ion component. The positive gradient exists at Ve < Vd 
and will amplify slow plasma modes with phase speed along the magnetic field w/fey < Vd- The 
dispersion relation of ion-sound waves gives u/k\\ = Cg (1 + fc^A^g)"^/^ and they are easily amplified 
when Vd> Cg. The excited ion-sound turbulence creates the anomalous resistivity and heating. 

With an increasing level of turbulence, C/turb/^^e^sTe, the mean free path of the electrons 
I decreases according to equation (53). This implies a higher resistivity and a lower thermal 
conductivity. Thus, a reduced i leads to a strong anomalous heating and the suppression of heat 
conduction out of the layer. It is then possible that the layer is overheated, loses hydrostatic 
balance, and expands in a runway manner. We now estimate the critical (. for this to happen. We 
neglect the radiative losses in this estimate. 

The ohmic heating rate per unit area of the layer is 

i^Ohm^/l-, (61) 

a 

where h = k-QT(./mig, a = nf.e^l/rrif.Ve is electrical conductivity, and = {k-BTe/nie)^^^ is electron 
thermal speed. The conductive heat loss Fc is given by equation (54), and we find 
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Hence the layer must be overheated if 

-^<2xl0-3 6'-3/2rc. (63) 

^Coul 

Here Tc = crTRj^ssj/ce ~ 10~^ — 10~^ (§ 6.2), and the temperature is limited by pair creation 
to ~ 1 (see § 6.3). One can see that the ohmic overheating requires a high level of plasma 
turbulence, when the mean-free path is strongly reduced compared to the Coulomb value. This 
level is achievable in principle because it is still below the upper limit C^turb/^e^B^^e ^ 1 that 
corresponds to the minimum i ^ Adb- Indeed, 

^^^,r^MV''flY'''r^(^_r,^y'\,, 10- 4^!. (64) 
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Note that anomalous heating becomes interesting only when the thermal heat flux is strongly 
suppressed. In other words, Fohm may compete with Fc and overheat the layer only when Fc is 
small. One can derive from the above formulae, 

rji / p \ —2 

Ohm / \ /rr\ 

(65) 



Fc \4 X 10^2 TcOe erg s~^ cm~2 
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Substituting Tc 



10-1 - 



10 ^ and ~ 1 one finds that the maximum energy flux from the 



ohmically-heated layer is ~ 10^^ erg s"^ cm~^. Even if the layer covers the whole star with 
area 47ri?^g ~ 10^^ cm^, its total power would not exceed 10'^'^ erg s^-*^. This is two orders of 
magnitude smaller than the observed luminosity of the corona, and hence the anomalous heating 
of the atmosphere is unlikely to generate the observed activity. 

The marginally stable state -Fohm = Fc has an electrostatic potential drop eS^^e ~ feB^e across 
the layer, which equals the energy gained by an electron crossing the layer. Since Tg is limited by 
pair creation, the maximum potential drop that can be sustained by the layer is ~ mgCp- This 
is at least 2 orders of magnitude lower than the voltage 7res'TieC^ that is developed in the corona 
with discharges. If the turbulence in the transition layer reaches the level required for runaway 
heating, the plasma may be supplied to the corona at a much lower voltage compared with the pair 
discharge described in § 5. 

Pair creation in the atmosphere (see § 6.3) would tend to reduce its mean molecular weight 
and increase h. When the expressions for Fohm and Fc are modified to include the pair density, 
ne — > rig + n±, one finds that the critical value of ^/^coui for runaway heating (eq. 63) increases as 
1 + n±/ne- The critical value of pair density is ^ nc{h/£), at which the diffusive fiux of hot 6=*= 
out of the atmosphere dominates the corona and shortens out the electric circuit. One finds the 
relation 



If the atmosphere is the main source of charges for the coronal current then n± ^ n* is expected. 
Ohmic heating is unimportant in such an "electrically regulated" atmosphere if n± !^ < Ue- It 
can become important if the atmosphere happens to enter the pair-dominated regime, n± > Ug- 
Then the overheating with -Fohm/^c (^e^ > 1 is possible, leading to over-production of pairs, 
n± > n*, which will suppress -fohm/-^c and quickly shut the instability. Perhaps a cyclic overheating 
with pair enrichment would occur. 

The effects discussed in this section indicate the possibility of an alternative way of forming 
the corona, through heating and pair creation in the atmosphere. This would require a high level of 
turbulence in the atmospheric plasma, and it is unclear whether this level is reached. Observations 
of the high-energy tail of the X-ray spectrum may help constrain the composition and temperature 
of the transition layer. 
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7. 



CRUST EXCAVATION AND FORMATION OF A LIGHT-ELEMENT 

SURFACE LAYER 



7.1. 



Mass Transfer in the Coronal Circuit 



The existence of an ion current in the magnetosphere implies the transfer of mass from the 
anode footpoint of a magnetic line to its cathode footpoint. The rate of transfer may be easily 
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estimated for the idealized configuration of a twisted dipole (§ 2), 

M - 47ri?^g ~ (^j^ ^ Sj BR^s ~ 10^^ (^j^ ejBi^ g s-\ (67) 

Here the dimensionless coefficient Sj = Acf) sin^ 9 measures the strength of the twist in the dipole 
magnetosphere (see eq. [1]). During the active life of a magnetar, t ~ 10'^ yr, a significant mass is 
transferred through the magnetosphere, AM 3 x 10^^ {ji/j) -Sis g- It is much larger than the 
mass of the stellar atmosphere, and the atmosphere is quickly removed from the anode part of the 
stellar surface. 

There is, however, a mechanism which can re-generate the atmosphere and regulate its column 
density. A thin atmosphere is unable to stop the energetic particles flowing from the corona and 
these particles will hit the crust and knock out new ions. In a steady state, the atmosphere column 
density can adjust so that the rate of ion supply equals the rate of their transfer through the corona. 
Then the rate of crust excavation at the anode footpoints is proportional to the ion current through 
the corona. 

Excavation depends on the chemical composition of the uppermost crust of the neutron star. If 
it is made of light elements, the atmosphere is easily re-generated by simple knock out of ions from 
the surface. If it is made of carbon, oxygen, or any heavier elements, the ions are condensed into 
long molecular chains in a 10^"'^ G magnetic field (Neuhauser, Koonin, & Langanke 1987; Lieb, 
Solovej, & Yngvason 1992; see Thompson et al. 2000 for estimates of the chain binding energy in 
the case of ~ 10^^ G magnetic fields) . The knock out of light elements by the beam is still possible 
in this case because of spallation reactions. We discuss the excavation process in more detail below. 



7.2. Excavation of Hydrogen and Helium Layers 

Light elements — hydrogen and helium — may be present only above a certain depth that 
may be estimated as follows. The column mass density E above a given depth in the crust is related 
to pressure P at this depth through the equation of hydrostatic equilibrium, 

ztt/?/ ztt 

Here it is taken into account that the pressure is dominated by degenerate electrons that are 
confined to their lowest Landau state {Ep < Eb in the density range of interest); pp ^ nieCp-jB^^ 
and Ep = \p\c^ + mlc^Y^'^ is the Fermi energy. A conservative upper bound on the depth of a 
hydrogen layer is obtained by balancing Ep with the threshold for electron captures on protons, 
Ep ~ 1.3 MeV. This corresponds to p ~ 10^ g cm~^. At larger depths, free protons are rapidly 
consumed by reaction e" + p ^ n + v. 

A more realistic upper bound takes into account that hydrogen is burned into helium, and the 
maximum depth of a helium layer is limited by the triple-alpha reaction 3a — >^^C, which is very 
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temperature-dependcnt. The interior of a magnetar probably sustains a temperature of ~ 3 x 10^ K 
at an age of ~ 10^ — IC^ yr (e.g. Arras, Gumming, & Thompson 2004). In the presence of a light- 
element atmosphere, the crust temperature remains close to the central temperature at depths 
p > 10^ g cm~^ (e.g. Potekhin et al. 2003). At such densities and temperatures helium (and 
hydrogen) cannot be present — its lifetime would be only ~ 10~^ (yePe)^ jr where Yg 0.5 is 
the electron fraction. Helium may survive only in an upper layer of density p < 10^ g cm~^, where 
temperature is below ~ 1 x 10^ K. 

The maximal layer of light elements will be drained rapidly from the anode surface and de- 
posited at the cathode surface in the presence of a magnetic twist. The excavation time down to a 
depth with a Fermi momentum pp is 

Z/A fji\~^f Pf Y -^NS,6 



ArUpji/e £j \j J \mecj gu 

Substituting ji/j ~ 0.1 and ej ~ 0.1 shows that the excavation time of the maximal hydrogen layer 
[Ep 2± 1.3 MeV) is much shorter than the active ~ 10^-yr lifetime of a magnetar. Taking a more 
realistic maximum density p ~ lO""* g cm~^ and pp/meC ~ 6 x pr,B^^ , one finds the excavation 
time of hours or days. It is faster than hydrogen consumption through diflFusion to a deeper layer 
of carbon (Chang, Arras, & Bildsten 2004). 

We conclude that the anode surface of the star should be made of heavier elements — carbon 
or iron. A gaseous atmosphere of light elements may be maintained above such a condensed, 
heavy-element surface due to the beam of energetic particles from the corona. 



7.3. Beam Stopping and Spallation of Nuclei in a Solid Crust 

When the coronal beam enters the crust, it is stopped by two-body collisions, and collective 
processes become negligible. In particular, the growth length of beam instability (eq. [46]) 
is much larger than the mean free path to two-body collisions (their ratio scales as v^J^ and is 
large in condensed matter because rie is large). The stopping length of the beam is then mainly 
determined by Coulomb collisions with the nuclei. 

Consider a downward-moving electron with an energy Ef, = ^rrieCp ^ Eb (eq. [29]). Collision 
with an ion can excite the electron to any Landau level n > with energy En = {2nB / Bqeb + 
1)"*^/^ < Ee- Most frequently, it will be the first level Ei = Eb because this process has the largest 
cross section (eq. 57). The probability of excitation to levels n' > n is reduced by a factor ~ ^ (it 
requires a smaller impact parameter; see also Kotov & Kelner 1985). One may assume a two-level 
system n = 0, 1 in simplified estimates. 

The excited electron moves along the magnetic line with a new Lorentz factor 7||, which is 
found from the energy conservation, 7||-E'b = Eg. It immediately de-excites and emits a photon at 
pitch angle ^em with energy Ej{9em) (§ 5.1). The normal modes of radiation in the uppermost crust 
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are determined by both vaccuum and matter polarization, which gives eUiptical polarization states. 
A photon in cither elliptical mode immediately converts to off the magnetic field if it is above 
the threshold, E-y > 2m(,c^ / sin^cm- This condition is satisfied for a fraction fs of de-excitation 
photons, which is estimated using the results of Herold et al. (1982). When B ^ i?QED, the 
angular distribution of the emitted photons is proportional to E"^ 

in both polarization modes. Combining this expression with Ej{6em) (eq. [36]), we find that fs 
can be well approximated by 

~ 1 - (-^) ~ 1 - 0-14 B^,\ B » i?QED. (70) 



(e - 1) V^QED, 

The net result of collisional excitation, de-excitation, and photon conversion to is the sharing of 
initial electron energy Ef. between three particles, and each of them will again collide with an ion 
and produce more e^. Thus, a pair cascade develops as a beam with initial 7 = 70 propagates into 
the crust, giving more with lower Lorentz factors 7 < 70- 

This branch of the cascade operates for a large fraction of the de-excitation photons. Photons 
that cannot immediately convert off the field (a fraction 1 — Jb) will wait until they interact with 
an ion in the crust (photon splitting is slow compared with photon- ion interaction) . Normally, such 
a photon will convert to in the Coulomb field of the nucleus. Less frequently, it can knock out 
a neutron, proton, or (rarely) alpha particle. This second branch of the cascade is most important 
for us, since it produces light elements. 

The probability of proton knockout instead of normal conversion to may be roughly esti- 
mated using the known nonmagnetic cross sections for photon-nucleous collision. Knockout requires 
a minimum photon energy of ~ 15 — 20 McV, and its cross section typically peaks at 20 — 30 
MeV (e.g. Dietrich & Berman 1988). The cross section for proton knockout from oxygen is 
(7(7, -"^^ O p,^^^ N) w 12 mb at Ej ~ 25 MeV. It scales approximately in proportion to Z, and so 
we take 

a{j,Z ^p,Z -1) = 1.5Z mb. (71) 
The cross section for pair creation is 



In 



2E^ \ 109 



rrieC^ J 42 



C7T « 3.2 Z"^ mb. (72) 



We can now estimate the number of protons knocked out per incident energetic electron. 
Wc will assume in this estimate an idealized surface layer that is composed of ions of a single 
charge Z. We roughly estimate the number of generations in the cascade down to an energy of 
~ 25 McV ~ 50meC^ as Ng ^ log2(7o/50) ^ 5. The number of ^ 25 McV photons created per 
incoming electron is about 7o/50. A fraction 1 — of these photons hit a nucleus, yielding a 
proton with a probability a (7, Z ^ p, Z — 1) /a^+^-. The resulting yield of protons per incident 
primary electron is 

. , _ (2^) -'^-^-"-^-^ ^ 0. (i) (I) . (73) 
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The cascade develops at a significant depth that buries the created positrons — they 
annihilate in the crust and cannot escape back through the surface. The first Coulomb collision 
occurs at a Thomson depth tt ^ Zo^ j (Jqo^x ~ I^Bi^^jZ. The relativistic cascade so initiated is 
strongly beamed into the star, and noticeable backscattering may occur only after the last stage 
of the cascade that gives with 7 ~ EB/'2,meC^ ~ "^^i^ 1 after log2(7omeC^/£^s) ~ 8 generations. 
Even the first branch of the cascade, which involves immediate pair creation off the magnetic 
field, ends at a significant tt ^ 100Z^^I?i5. This layer is optically thick to positron annihilation. 
Therefore, the anode surface can supply positive charges only by emitting light ions, mostly protons 
that have been produced by spallation. 

Spallation reactions occur down to a Thomson depth tt ~ 
knocked-out protons diffuse through the spallation layer on a short timescale, much shorter than 
the excavation time of this layer, and are quickly sublimated to the atmosphere. Thus, excavation 
of the crust is driven by spallation followed immediately by sublimation. 

Our preferred circuit model has voltage e$e ~ 1 GeV and ion current jj ~ O.lj (§ 5). If the 
gaseous atmosphere does not decelerate the particles flowing from the corona, they will impact the 
solid crust with 70 ~ 10^. The estimated proton yield Np ~ 0.1(.Z/26)~^ is then sufficient to feed 
ji ~ O.lj and can supply more ions than transferred through the corona. Then the column density 
of the atmosphere will grow until it is able to damp the relativistic beam from the corona and 
reduce the crustal excavation rate. 

The profile of the atmosphere is determined by its efficiency of stopping the relativistic e^. 
If collisionless processes do not decelerate the relativistic beam but just form a plateau in its 
distribution function (§ 6.1), then a thick atmosphere must be accumulated to decelerate the beam 
through Coulomb collisions. Collisions initiate a cascade similar to that in the crust: the beam 
particles collide with ions and emit photons that convert to pairs. The crustal spallation rate will 
be reduced if the beam is decelerated to 7 ;^ 50 before reaching the solid surface. This requires 
log2(7res/50) ~ 5 generations of pair-creation, and so the optical depth to Coulomb collisions should 
be ~ 5. Thus, the column density of the atmosphere will be regulated to the value 

log2(7rcs/50) 1 , . 

J^atm ~ ^ lUU -Di5 (7^ . (74) 

f^Coul 

If Np/{1 + Np) < ji/ j, the crust will not supply the ion flux implied by our circuit solution 
with e$e ~ lrcsfn.e(?- Wc then expect the voltage to grow, giving a higher spallation rate of the 
crust. Wo leave this regime for a future study. 



^"Protons have a binding energy Ei of a few keV and are easily sublimated at the temperature k^T ~ 0.5 keV, in 
contrast to the heavy elements that are locked in molecular chains with Ei ~ 10^ keV. 

^^The relatively low density of the atmosphere implies two differences with the cascade in the crust: (1) The 
polarization states of photons are now linear — vaccuum O- and E-modes - and the E-mode does not convert to 
off the field. (2) Splitting of E-mode occurs faster than photon-ion interaction. This leaves only one branch for the 
cascade, with production off the magnetic field. 
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The thickness of the plasma layer that is accumulated near the cathode surface is limited by 
burning processes and diffusion to greater depths (Chang et al. 2004). 

8. DISCUSSION 

The basic finding of this paper is that an corona must be maintained around a magnctar. 
It exists in a state of self-organized criticality, near the threshold for breakdown. The stochastic 
discharge continually replenishes the coronal plasma, which is lost on the dynamic time tdyn ~ 
10-4 s. 

The solution for the plasma dynamics in the corona is strongly non- linear and time-dependent. 
It is essentially global in the sense that the plasma behavior near one footpoint of a magnetospheric 
field line is coupled to the behavior near the other footpoint. Remarkably, this complicated global 
behavior may be described as essentially one-dimensional electric circuit that is subject to simple 
boundary conditions E ^ on the surface of the star and can be studied using a direct numerical 
experiment. 

The established voltage across the entire magnetospheric circuit is marginally sufficient to 
accelerate an electron (or positron) to the energy where collisions with ambient X-ray photons 
spawn new pairs. This physical condition determines the rate of energy dissipation in the corona. 
The voltage is maintained by small deviations from charge neutrality in the coronal plasma as 
required by Gauss' law. The generation of this voltage is precisely the self-induction effect of the 
gradually decaying magnetic twist, and the energy release in the corona is fed by the magnetic 
energy of the twisted field. 

8.1. Comparison With Canonical Radio Pulsars 

Canonical radio pulsars are usually assumed to have a magnetosphere that is almost everywhere 
potential, V x B = 0, except for narrow bundles of open field lines that connect the polar caps of the 
star to its light cylinder. The plasma content of the magnetosphere then depends on its rotation: 
the charge density pco = ■ B/27rc must be maintained to screen the electric field induced by 
rotation (Goldreich &; Julian 1969). In this picture, the persistent current of a canonical pulsar 
only flows along the open field lines out to the light cylinder. The sign of this current /open is 
determined by the corotation charge density. The closure of the electric circuit back to the star 
remains poorly understood. 

Previous studies of plasma near neutron stars focused on open magnetic lines. It was shown 
that pair discharge can be sparked near the polar cap (Rudcrman & Sutherland 1975; Arons & 
Scharlemann 1979; Muslimov &: Tsygan 1992), and the created pairs carry the current /open toward 
the light cylinder. Recent calculations of the discharge (Harding & Muslimov 2002) employ an 
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approximate analytic prescription for the decay of the electric field with distance from the star. 
(Self-consistent calculations of the electric field and pair creation are now available for the outer 
gap, sec Takata ct al. 2006.) The global dynamics of the plasma on open magnetic lines is 
poorly understood because of the complicated behavior at the light cylinder, so only one boundary 
condition is known — at the polar cap. 

The solution presented in this paper applies to currents along closed magnetic lines of the inner 
magnetosphere. In this case, it is easier to find the global solution for the plasma dynamics because 
both boundary conditions are defined at the two footpoints of the magnetic lines. This problem 
has important applications to the magnetar physics for two reasons: (1) The observed radiative 
output is generated in the closed magnetosphere — the dissipation rate on the open field lines is 
orders of magnitude smaller and cannot explain the observed emission. (2) A current is naturally 
created on closed magnetic lines because they are twisted by the starquakes of magnetars, and our 
solution shows that such twists have a long lifetime. 

One may compare the twist of open magnetic lines at the polar cap, as sustained by the rotation 
of the star, B^/Bp ~ (i?Ns/-Ric)^! with the twist that can be created in the closed magnetosphere 
by the starquakes. Consider a current / flowing along a magnetic flux tube of radius a and let B^ 
be the twisted component of the field. Then 

c c o do!) , , 

I = -^aB, = -a-Bp^^, (75) 

where I is the length measured along the tube. Since the current / and the poloidal magnetic flux 
na^Bp are constants along the tube, the twist angle satisfies d^/dZ = const. Hence, 

h7 ~ ~R — ' 

where Rmax is the maximum radius to which the tube extends and Acj) is the net twist angle 
integrated along the tube. The net twist angle is limited by the global stability of magnetosphere 
to A(p ^ 1. On the star's surface one then finds B^/Bp = ad(j)/dl ^ a/R^g^x- This implies that 
the closed field lines with -Rmax ~ 2i?NS and a ~ i?NS can bear much stronger twists than the open 
lines. 

Note also one more differrence of the magnetar currents: they impact the stellar surface with a 
modest Lorentz factor 7^ ~ 10^. Such a beam can deposit its energy in an optically thin atmospheric 
layer through a plasma instability (§ 6.1). By contrast, from the discharge in an ordinary radio 
pulsar impact the star with high 7^ = 10^ — 10^, which suppresses the beam instability in the 
atmosphere: its length scale li, is proportional to 7;, (eq. 46) and becomes much larger than the 
atmosphere scale-height. Such a beam must be stopped by Coulomb collisions and relativistic 
bremsstrahlung in the solid crust. Therefore, the polar cap of an ordinary pulsar is expected to 
release the beam energy in the form of blackbody radiation, in agreement with observations (e.g. De 
Luca et al. 2005). By contrast, magnetars may have a strongly overheated hydrostatic atmosphere, 
with kBT > 100 keV (§ 6.2). 
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8.2. Comparison with Solar Corona 

The corona of a magnetar may be viewed as a collection of closed current-carrying flux tubes 
with footpoints anchored in a relatively cold stellar surface. As in the Solar corona, a current is 
created by the twisting motions of the footpoints. Another similarity with the Solar corona is the 
important interaction of the coronal plasma with a dense atmospheric layer on the stellar surface 
where the plasma becomes coUisional and can efficiently radiate the energy flux received from the 
corona. 

However, the magnetar corona is distinguished from the Solar corona in a few important 
respects. The differences in the behavior of the magnetic field are as follows. 

1. The footpoints of magnetic lines of the Solar corona are constantly pushed and twisted because 
the outer envelope of the Sun is convective and in constant motion. In addition, magnetic flux 
tubes are observed to emerge from the Solar convection zone with a twist already implanted in 
them. By contrast, the neutron star's rigid crust is subject only to rare and sporadic yielding 
events (starquakes) . Because the crust is stably stratified, the transfer of a toroidal magnetic field 
to the exterior occurs through a differential rotation of the crustal material along the gravitational 
equipotential surfaces. 

2. The magnetic twists have long lifetimes between starquakes, and their evolution differs from 
the evolution of twists in the Solar corona. Once a current j = (c/47r)V x B is created in the 
magnetospherc of a neutron star, it lives a long time because the net current / ~ -Rns.? enormous 
and maintained by its self-induction. The gradual decay of the current induces an e.m.f. along the 
magnetic field lines, which forces the plasma and electric field to self-organize into a quasi-steady 
electric circuit. The ease with which charges can be supplied by 6=*= production when the induced 
voltage exceeds e$e ~ 1 GeV implies a "bottleneck" for the decay of the twist (§§ 2 and 5). 

3. The evolution of the magnetic field between the starquakes involves a gradual release of the twist 
energy through Joule dissipation and spreading of the twist toward the magnetic polar axis. The 
latter process causes a gradual flaring of the poloidal field lines. Explosive magnetic reconnection 
similar to Solar flares is likely to be triggered by starquakes, but not during the following slow- 
dissipation stage. 

4. The external magnetic field of a magnetar is likely to have less small-scale structure compared 
with the Solar corona. The small-scale braiding of magnetic lines followed by reconnection and 
dissipation is one of the possible heating mechanisms for the Solar corona (Parker 1979). Precisely 
what creates this structure is unclear: it could be driven by the hydrodynamic turbulence that is 
excited at the top of the Solar convection zone or by the Taylor relaxation of a helical magnetic 
field (e.g. Diamond & Malkov 2003). In magnetars, the small-scale structure of magnetic field can 
be expected to dissipate quickly after the starquakes if it is created in these events. The dissipation 
timcscalc is longest for the large-scale twists expected from large-scale motions of the crust. Indirect 
evidence that small-scale instabilities are less important comes from the observed distribution of 
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energies of SGR bursts: the cumulative energy in bursts E^dN/dE ~ E^-^ is weighted toward the 

largest events (Cheng et al. 1996; Gogus et al. 2001). By contrast, several recent measurements 
of the distribution of Solar flare energies suggest a power- law dN/dE = E^ with k < -2, and 
are therefore consistent with the basic picture of coronal heating by 'nano- flares' (Walsh & Ireland 
2003). 

Let us also note a few differences in the plasma behavior. 

5. The rapid cooling of the transverse motion of charged particles in the ultra-strong magnetic field 
of magnetars makes the particle distribution functions one-dimensional. This reduces the number 
of relevant plasma modes in the corona. 

6. The coronal plasma of magnetars is continually lost to the surface and re-generated. By contrast, 
the Solar corona is in hydrostatic equilibrium. Its density n is much higher than the minimum 
needed to maintain the current, and the mean drift speed j/en of the coronal currents is much 
below the electron and ion thermal speeds. So, the mean current does not experience a strong 
anomalous resistivity. Strong plasma turbulence is expected in localized current sheets and can 
develop explosively, giving rise to the Solar flares. In the magnetar corona, a strong turbulence is 
maintained routinely, everywhere where current flows. The very mechanism of the current supply 
involves a strongly non-linear phenomenon - discharge. 

7. The magnetar corona is dominated by relativistic pairs. As indicated by our numerical 
experiments, ions lifted from the surface contribute a modest fraction, perhaps ~ 10%, to the 
charge flow. 

8.3. Magnetospheric Circuit as a Double-layer Problem 

Prom the plasma physics point of view, there is an interesting connection between the mag- 
netospheric circuit and the double-layer problem that was first studied by Langmuir (1929) and 
later by many authors (for reviews see Radau 1989; Block 1978). The magnetospheric current flows 
along a magnetic tube with footpoints on a cold conductor and is subject to the boundary condition 
E' ~ at the footpoints. We argued that the tube curvature is not essential to understanding many 
features of the circuit, and studied the charge flow in a straight tube where gravity is replaced by a 
force mg that attracts plasma to the anode and cathode boundaries. The longitudinal electric field 
must develop between the two boundaries because of gravity. This is a significant difference from 
the normally invoked mechanisms for double-layer formation (such as current-driven instability). 
In our case, gravity traps the thermal particles at the magnetic footpoints and impedes the flow of 
the current that is required by the self-induction of the twisted magnetic field. 

Thus, perfect conditions for double-layer formation are created. And indeed, our numerical 
experiment showed the formation of a huge relativistic double layer between the footpoints, of size 
L ~ -Rns- This result is obtained if the current can be maintained only by particles lifted from 
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the anode and cathode surfaces. It could have dramatic consequences for the magnetar physics 
because of the huge dissipation rate in the double layer: the energy of magnetic twist would be 
quickly released in a powerful burst. 

However, there is one more ingredient in the problem, which again makes the magnetar circuit 
different from the previously studied double layers: new charges (e^) can be created by accelerated 
particles in the tube. This process qualitatively changes the problem because voltage is now limited 
by discharge (§ 5). As a result, the circuit enters the new state of self-organized criticality where 
voltage repeatedly builds up and is screened by the creation. This allows the current to flow 
and maintain the required V x B. 

8.4. Observational Implications 

Two important observational implications of the magnetospheric twist were already pointed 
out by TLK, based on the force- free model of a twisted dipole: 

1. The poloidal field lines are inflated compared with a pure dipole, so that a larger magnetic 
flux connects the star to its light cylinder Ri^. This implies a higher spin-down torque. Thus, the 
twisting leads to a higher spindown rate P. 

2. The minimum particle density needed to maintain the currents, ric = j/ec, is high enough to 
make the magnctosphere semi-opaque. The main source of opacity is resonant cyclotron scattering. 
If the magnetic field lines arc twisted, then the optical depth to scattering happens to be ~ n/ric 
(independent of r, B, and nature of the flowing charges). Scattering by the magnetospheric charges 
must modify the X-ray spectrum and pulse profile. 

TLK also pointed out that if the magnetospheric current is carried by ions and electrons lifted 
from the stellar surface by a minimum voltage e$e ^ gmpRNS > then one expects a dissipation rate 
comparable to the observed persistent luminosity of magnetars. 

Further implications of our model for the coronal electric circuit are as follows. 

3. The density of the plasma corona is found to be close to its minimum ric = j/ec as conjectured 
by TLK. It is largely made of pairs rather than electron-ion plasma. This suggests that the 
optical depth of the corona is determined mostly by the cyclotron resonance of 6=*=. All species of 
particles {e~ , e+, and even ions) flow with relativistic bulk speeds and large temperatures ~ e$e- 
The cyclotron resonance is very broad, which will have important implications for radiative transfer 
through the corona. Note, however, that our calculations focused on small radii r ~ i?NS where 
most of the energy release occurs. At larger radii r ~ (10 — 20)i?NS) where electron cyclotron energy 
is in the keV range, feel a strong radiative drag force and can decelerate to a small velocity. We 
defer a detailed analysis of this part of the corona to another paper. 

4. The voltage e$e — 7res"^eC^ ^ 1 GeV along the twisted magnetic lines of the corona is regulated 
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by the threshold for discharge, and its exact value depends on B and the spectrum of target 
photons involved in creation, especially on its high-energy tail at 10 — 300 keV (§ 5). By 
coincidence, the mean electric field developed along the magnetic lines, E ~ ^^/R^Si is marginally 
sufficient to lift ions from the anode footpoint. Therefore, in addition to pairs, some ions are 
present in the corona. Their abundance depends on the exact $e ^'^d is about 10% for e$e = 1 GeV 
(see Fig. 10). The coronal ions may produce interesting features in the magnetar spectrum and 
likely to play a key role in the low- frequency emission (§ 8.4.2). 

5. The rate of energy dissipation in the twisted magnetosphere is given by L^iss = I^e where / is 
the net current through the corona. The current may be estimated as 

I^jsa^-^^Acp-^a^, (77) 
47r ii^NS 

where a is the size of a twisted region on the stellar surface and = |V x B|(5/i?Ns)~^ ~ 1 
characterizes the strength of the twist. The calculated e$e is comparable to 1 GeV and implies 

L„„ = ^ 10- A, (^) (^) ' (g^) erg s"'. (78) 

The observed luminosity i^diss ~ 10^^ erg s~^ is consistent with a partially twisted magnetosphere, 
a ~ 0.3i?NS; or a global moderate twist with A(^(e$e/GeV) ~ 0.1. (Larger dissipation rates and 
larger twists are not yet excluded by the data: e.g. the output at 1 MeV is presently unknown.) 
More precise formulae may be obtained for concrete magnetic configurations. In the simplest case of 
a twisted dipole (TLK), one can use equation (1) and find / by integrating ydS over one hemisphere 
of the star. This gives 1 = 1 A(f)cBpo\cRNS and Ldiss = 1-2 x lO^''A0Spoic,i5^NS,6^e,GeV- 

6. Once created, a magnetospheric twist has a relatively long but finite lifetime. The energy stored 
in it is ftwist ~ (-^^/c^)-^NS- This energy dissipates in a time 

/ \ / \ 2 



"""""^ ~ ^ ~ " c2 $2 ^ Vl036 erg s-1 j ) ^ ^ 

The exact numerical coefficient in this equation again depends on the concrete magnetic configu- 
ration. For example, the energy stored in the global twist of a dipole is ftwist = O.17(A0)2ii^(jjpoig 
(TLK), where -Edipoic = n ^poie-^NS energy of a normal dipole. In this case, one finds 

idecay = 0.8 (i?Ns/c^) ^diss/^e- 

The characteristic timescale of decay is determined by voltage and the corona luminosity 
-^^diss- Note that a stronger twist (brighter corona) lives longer. Given a measured luminosity and 
its decay time, one may infer $e and compare it with the theoretically expected value. Since ^>e 
does not change as the twist decays in our model, one can solve equation (79) for £^twist(i) and 
-^diss(i) = -d^twist/dt. This gives 

LdiUt) = 5 X 1035 (-^^ {to - t) erg s'S (80) 
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where time t is expressed in years. The final decay of the twist should occur at a well-defined time 

to = t + tdccay • The finite time of decay implies that the emission called "persistent" in this paper 
is not, strictly speaking, persistent. If the time between large-scale starquakes is longer than tdccayj 
the magnetar should be seen to enter a quiescent state and the observed luminosity should then be 
dominated by the surface blackbody emission. 

7. The voltage <&e implies a certain effective resistivity of the corona, TZ = ^e/I- This resistivity 
leads to spreading of the electric current across the magnetic lines. Equivalently, the magnetic 
helicity is redistributed within the magnetosphere. This spreading is described by the induction 
equation d'B/dt = — cV x E. The magnitude of V x E depends on the gradient of -Ey transverse 
to the twisted magnetic field lines. Our model of the magnetospheric voltage implies that $e has a 
gradient on a lengthscale ~ i?NS) so that dBtwist/dt ~ cE^^/R^s. The characteristic spreading time 
of the twist -Btwist = A(j)B (a/i?Ns) is ^spread ~ -Btwist [9 Bx^ist / dt)~^ and can be roughly estimated 
as 

If a twisting event happens near the magnetic equator, it initially does not affect the extended 
magnetic lines with i?max ^ R\c- Its impact on spindown may appear with a delay as long as 
^spread- Although the exact numerical coefficient in equation (81) is unknown, the rough estimate 
suggests a timescale that is comparable to the decay time (eq. [79]). Their ratio is oc (a/i?Ns)~^ 
and depends sensitively on geometrical factors. It is possible that starquakes occurring far from 
the poles will have little effect on the spindown. 

8. Even a small ion current implies a large transfer of mass through the magnetosphere over the 
^ 10^-yr active lifetime of a magnetar. Essentially all hydrogen and helium that could be stored 
in the upper crust is rapidly transferred from the anode surface to the cathode surface, thereby 
exposing heavier elements (e.g. carbon) that are strongly bound in molecular chains. A continuing 
flux of protons can, nonetheless, be maintained in the circuit. The bombarding relativistic electrons 
from the corona cause the gradual disruption of ions just below the solid surface. They trigger a 
cascade of e"^ pairs and gamma-rays that knock out nucleons and occasionally a-particles from the 
heavier nuclei. This process re-generates the light-element atmosphere on the surface and regulates 
its column density (§7). As a result, the Thomson depth of the atmosphere equilibrates to a 
value ~ 100. The chemical composition of the upper crust (carbon or heavier elements) and its 
light-element atmosphere affect radiative transfer and the emerging spectrum of the star's radiation 
(e.g. Zavlin & Pavlov 2002). 

9. Our study of the transition layer between the corona and the star (§ 6) suggests that the energy 
dissipated in the corona is radiated below it in two forms (see also Thompson Sz Beloborodov 
2005). First, the high-energy radiation with a hard spectrum possibly extending up to ^ 1 MeV; 
this radiation is caused by collisionless dissipation of the coronal beam in the transition layer. 
Second, a blackbody radiation component that is caused by thermalization of the remaining energy 
of the beam when it enters the dense atmosphere and then the crust; two-body collisions and 
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development of an e cascade are responsible for this thermalization. A crude estimate gives 
comparable luminosities in the high-energy and blackbody components. Further investigation of 
the transition layer may show that the collisionless dissipation is more efficient. This would be in a 
better agreement with observations of SGR 1806-20 whose high-energy component was found more 
luminous than the blackbody component (Mereghetti et al. 2005; Molkov et al. 2005). 

8.5. Emission from the Magnetar Corona 

8.5.1. Hard X-ray Emission 

The radiative output from the corona is observed to peak above 100 keV (Kuiper, Hermsen 
& Mendez 2004; Mereghetti et al. 2005; Molkov et al. 2005; den Hartog et al. 2006). It is likely 
fed by the inner parts of the corona: in the model investigated in this paper, one expects that the 
bulk of the coronal current is concentrated fairly close to the star, at radii < 2i?NS- Currents in 
the outer region r :§> i?NS (even if it is strongly twisted) flow through a small polar cap on the star 
and therefore are small (§ 8.1). 

Emission from the inner corona must be suppressed above ^ 1 MeV, regardless the mechanism 
of emission, because photons with energy ^ 1 MeV cannot escape the ultra-strong magnetic field. 
There are two possible linear polarizations of photons in the corona: E-va.o(ie (electric vector 
perpendicular to B) and 0-mode (electric vector has a component along B). An energetic E- 
mode photon cannot escape because it splits into two photons. An 0-mode photon propagating 
at angle with respect to B cannot escape when its energy E^ > 2meC^/ sin^ because it converts 
to a pair. Therefore, the spectrum escaping the inner corona is limited to ~ 1 MeV. The available 
data show that the high-energy radiation extends above 100 keV. There is an indication for a cutoff 
between 200 keV and 1 MeV from COMPTEL upper limits for AXP 4U 0142-^61 (den Hartog et 
al. 2006). 

The possible mechanisms of emission are strongly constrained. The corona has a low density 
n ^ Uc = Jb/^c ^ B^Air eRjsfs)^^ ~ lO^^Sis cm^"^ and particle collisions are rare, so two-body 
radiative processes are negligible. The only available emission process in the corona is upscattering 
of surface X-rays. The scattering can occur resonantly (excitation to the first Landau level followed 
by de-excitation) or non-resonantly (Compton scattering). The scattering, however, cannot explain 
the observed 100-keV emission: 

1. Resonant upscattering by electrons (or positrons) with Lorentz factor 7e gives high-energy pho- 
tons with energy ~ je^eB/nieC » 1 MeV (§ 5.1), which cannot escape. They either split (if 
E-mode) or convert to pairs (if 0-mode). Splitting may give escaping photons below 1 MeV, but 
their number is not sufficient to explain the observed emission. The voltage along the magnetic 
lines adjusts so that the coronal particles experience ~ 1 resonant scattering before they are lost 
to the surface, while the observed emission requires the production of 10^ — 10^ 100-keV photons 
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per particle. 

2. Resonant upscattering by ions with Lorentz factors 7^ gives photons of energy 'yieB/nipC, which 
is too low to explain the 100-keV emission (ions are mildly relativistic in the circuit with ^ GeV 
voltage, 7i ~ 1). 

3. Non-resonant Compton scattering is negligible for £'-mode photons. The cross section of 0-mode 
scattering by electron with Lorentz factor 7e is ~ (JT/7e.^^ This may be seen by transforming the 
photon to the rest frame of the electron. Then the photon propagates at angle sin^' ~ l/7e with 
respect to the magnetic field and its electric vector is nearly perpendicular to B. Viewing it as a 
classical wave that shakes the ID electron, one finds the cross section of scattering a = ut sin^ 9' = 
(TT/7e- The corresponding optical depth of the corona is r = nc-RNSO" = Tc/le- The upscattered 
photon is boosted in energy by 7g and the net Compton amplification factor of 0-mode radiation 
passing through the corona is r7g = Tc, i.e. does not depend on je- Since Tc = cjt-Rns^c is small, 
~ 0.1 — 0.01, one concludes that the energy losses of the corona through non-resonant Compton 
scattering are a small fraction of the surface 0-mode luminosity. 

A possible source of observed 20 — 100 keV emission is the transition layer between the corona 
and the thermal photosphere (sec Thompson & Bcloborodov 2005 and § 6). This layer can radiate 
away (through bremsstrahlung) a significant fraction of the energy received from the corona. Its 
temperature is ~ 100(^/^coui)~^^^ keV (eq. [56]), where i/lcoui < 1 parameterizes the sup- 
pression of thermal conductivity by plasma turbulence. This suppression is likely to be significant 
and increase the temperature of the layer. The emitted hard X-ray spectrum may be approximated 
as single-temperature optically thin bremsstrahlung (Thompson &: Bcloborodov 2005). Its photon 
index below the exponential cutoff is close to —1, in agreement with magnetar spectra observed 
with INTEGRAL and RXTE. 

The observed pulsed fraction increases toward the high-energy end of the spectrum and ap- 
proaches 100% at 100 keV (Kuiper et al. 2004). If the hard X-rays are produced by one or two 
twisted spots on the stellar surface, the large pulsed fraction implies that the spots almost disappear 
during some phase of rotation. Then they should not be too far from each other and, in particular, 
cannot be antipodal because a large fraction of the stellar surface is visible to observer due to the 
gravitational bending of light: S'vis/47ri?^g = {2 - AGM/c^Rf^s)'^ 3/4 (Bcloborodov 2002). The 
decreasing pulsed fraction at smaller energies may be explained by contamination of the strongly 
pulsating component by another less pulsed and softer component of the magnetospheric emission 
(see § 8.4.2). 

Hard X-ray emission may also be produced at about 100 km above the star surface (Thompson 
Sz Beloborodov 2005). The magnetic field is about 10^^ G in this region, and the electrons can 
experience a significant radiative drag as a result of resonant cyclotron scattering. They cannot close 



^■^This assumes Thomson regime of seattering, jefitu <ti rne(? ■ At larger energies, the electron recoil becomes 
important (Klein-Nishina regime) and the cross section is reduced. 
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the circuit unless an additional electric field develops that helps them flow back to the star against 

the radiative pressure. Then some positrons in this region can undergo a runaway acceleration and 
Compton scatter keV photons to a characteristic energy Ery ~ {3-ir/8)a~^mcC^ , where a = is 
the fine structure constant. Such a photon is immediately converted to an pair in the magnetic 
field, and the created pairs emit a synchrotron spectrum consistent with the observed 20-100 keV 
emission. 

8.5.2. 2-10 keV Emission 

A promising emission mechanism for the observed soft X-ray tail of the blackbody component 
is resonant upscattering of keV photons in the corona (TLK). The emergent soft X-ray spectrum 
and pulse profile is likely to form at radii where B 10^^ G and the can be decelerated to 
mildly relativistic speeds by the resonance with keV photons. 

Ions flowing through the inner magnetosphere can also upscatter the keV photons. The ions are 
mildly relativistic and their cyclotron resonance near the star is 6i3i5 keV — just in the 2 — 10 keV 
band. On the other hand, the ion current may be suppressed (see Fig. 10), which would imply their 
small optical depth in the corona, except the first 10^ — 10^ seconds following an SGR flare when 
many ions are lifted to the magnetosphere (Ibrahim et al. 2001). 

8.5.3. Optical and Infrared Emission 

The observed infrared (K-band) and optical luminosities of magnetars are ~ 10^^ erg/s (Hulle- 
man et al. 2004), which is far above the Rayleigh- Jeans tail of the surface blackbody radiation. 
The inferred brightness temperature is ~ 10^'^ K if the radiation is emitted from the surface of the 
star. Four principle emission mechanisms may then be considered. 

1. Coherent plasma emission. Eichler, Gedalin, & Lyubarsky (2002) proposed that the optical 
emission is generated as 0-mode radiation at a frequency ~ wpg. We find that this mechanism 
gives too low a frequency of emission. The plasma frequency in the corona is suppressed by the large 
relativistic dispersion of the e^, tope ~ (47rnce^/me7e)^''^; in addition, it may not be blueshifted by 
bulk relativistic motion as suggested by pulsar models. A fraction of ~ photons may, however, 
be upscattered to the optical or higher bands by the coronal that have a broad distribution 
1< 7e < 10=^. 

The plasma frequency is higher in the hydrostatic atmosphere, where the density is higher. For 
a hydrostatic layer of temperature Tg and Thomson depth tt one finds fkope = 0.65 {TTrrieC^ / kT^)^^'^ eV. 
In this case, however, the brightness temperature of the 0-mode photons is limited to a value 
^B^fe ~ me,(? /tt by induced downscattering in the atmosphere. The corresponding upper limit for 
K-band luminosity, ujLi^ = (kBTbCV^ /Sn^c^) AttR"^^ ^ 10^^ Tr^^ ('^/'^K-band)'^ erg s~^, is in conflict 
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with observations. 

2. Non-thermal synchrotron emission from electrons with high 7e. This mechanism could work 
only far from the star, r ^ 200i?NS7e^^) where the minimum energy of the synchrotron photons, 
~ 'ygheB/nieC, is below the optical band {B < 10^7~^ G). The magnetic lines extending to large 
r form a small cap on the stellar surface and carry a tiny fraction of the total magnetospheric 
current, / ~ (''/-Rns)"^ ~ 3 x 10~^7e If the voltage maintained along the extended lines is 
comparable to 1 GeV, energy dissipation rate in the region is only 10^^77'^^^ erg s^^, orders of 
magnitude smaller than the observed optical-IR luminosity. Therefore this mechanism appears to 
be unlikely. 

3. Ion cyclotron emission from radii ~ 20i?NS) where the ion cyclotron frequency is in the optical 
band. The observed flux then implies ion temperatures of ~ 10^^ K. Such temperatures are indeed 
achieved in the corona. This mechanism also requires a sufficient number of ions in the corona, 
which is sensitive to the ratio of voltage e$e to the ion gravitational energy (Fig. 10). The ions 
will act as a filter for the radio and microwave radiation that is emitted by coherent processes 
near the star. The radiation reaching r > 500 km will be partially absorbed by the ions at their 
cyclotron resonance, pumping their perpendicular energy E±. This energy further increases as the 
ions approach the star since E±/B is an adiabatic invariant. Their cyclotron cooling and transit 
times become comparable at a radius r* ~ 200 (vj/c)"^/^ -^ns 6 -^pofe 15 Then the ions emit 
cyclotron radiation in the optical-IR range, since heB{r*)/mpC -Spoie^i5(i'i/c)^/^i?i^g^g eV. 

4. Curvature emission by pairs in the inner magnetosphere may extend into the optical-IR range. 
Charges moving on field lines with a curvature radius Rq ^ i?NS will emit radiation with frequency 

~ 7fc/27ri?c', which is in the K-band {u = 10^'^ Hz) or optical when -jemeC^ ^ GcV. As 
in models of pulsar radio emission, bunching of the radiating charges is required to produce the 
observed luminosity (e.g. Asseo, Pelletier, & Sol 1990). The radio luminosity of a normal pulsar 
can approach ~ 10~^ of its spindown power, and it is not implausible that the observed fraction 
~ 10~^ of the power dissipated in a magnetar corona is radiated in optical-IR photons by the same 
coherent mechanism. 

There are indications of charge bunching on the Debye scale in our numerical experiment, 
although this detail of the experiment is difficult to scale to the real corona of magnetars. Suppose a 
fraction / of the net flow of particles through the corona Nf. = I / e ^ lO'^^Ljissse (e$e/GeV)^^ s^^ 
is carried by bunches of M electrons (or positrons). The power emitted by one bunch is -B = 
(2/3)(A/'e)'^7gC/i?^, and the total coherent luminosity produced by the corona is given hy L k, 
fAfNeJe e^/Rc- Using the relation uc ~ (c/27ri?c') 7g for curvature radiation, we find the number 
of charges per bunch that is required to explain the observed luminosity L ~ 10^^ erg s"^, 

J\f ~ 10^ Ls2i^cM. (82) 

This is much smaller then the number of charges per Debye sphere, A/d ~ ('^c''e)~^^^ ~ lO"*"^. Thus, 
a very weak bunching of charges (or very small bunches, with a scale ~ 10~^A£)e) can generate a 
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significant power in coherent curvature radiation. The required bunching is small compared with 
what is commonly invoked to explain the emission of radio pulsars. 

For example, one can compare with the Crab pulsar. The total flux of pairs from the pulsar, 
iVg ~ 10^^ s~^, is comparable to the flux in a magnetar corona. The pulsar flow is concentrated 
on the open field lines, and its density is enhanced by a factor i?ic/^NS ~ 10^. The curvature 
radius of the open field lines near the star is Rc ~ 10i?NS- These differences are, however, modest 
compared with the much larger bunching efficiency that is implied by the observed brightness of 
the radio emission from the Crab pulsar. Supplying its luminosity L ~ 10^^ erg s~^ at ~ 10^ Hz 
would require 7e ~ 60 and M ~ 10^^/~^. 

The dispersion relation for the IR-optical radiation in the corona is dominated by plasma 
rather than by vacuum polarization. The curvature radiation is, generally, a superposition of all 
possible polarization modes (see Arons & Barnard 1986 for a description of the normal modes of 
a strongly magnetized plasma). Emission of the superluminal 0-mode will occur if the bunch size 
(measured in its rest frame) is smaller than the plasma skin depth. The E-mode also comprises a 
fraction of the curvature emission when the magnetic field lines are twisted (Luo & Melrose 1995). 
Curvature emission in the subluminal Alfven mode is also possible, but this mode is ducted along 
the magnetic field lines and has a frequency below the ambient plasma frequency. 

The spectrum of curvature emission is expected to cut off at frequencies above vq that cor- 
responds to the maximum Lorentz factor (a few 7res) of in the corona. This cutoff may be 

observable in the optical or infrared band (cf. Hulleman et al. 2004). Note also that the K-band 
brightness temperature inferred from observations of AXPs is remarkably close to the effective tem- 
perature of the coronal e^, T^q ^ 7emeC^/A;B, which is bounded by ~ 10^"^ K by the pair creation 
process. 

8.6. Evolution Following Bursts of Activity 

Bursts of magnetar activity are associated with starquakes — sudden crust deformations that 
release the internal magnetic helicity of the star and impart an additional twist to the magneto- 
sphere. The increase of the twist should lead to a larger luminosity from the corona since the 
luminosity is simply proportional to the net current flowing through the magnetosphere. Between 
bursts, the magnetospheric currents persistently dissipate and the corona disappears after the time 
idecay (eqs. [79] and [80]). If there were no bursts during a long time t > idecay t^ie magnetar 
should enter the quiescent state with no coronal emission. Such behavior has been observed in 
AXP J1810-197 (Gotthelf & Halpern 2005). An outburst, which occured between November 2002 
and January 2003, was followed by the gradual decay of soft X-ray emission on a time-scale of a 
few years, and now the source is returning to the quiesent state. It would be especially interesting 
to observe the evolution of the hard X-ray component, which is not detected in the relatively weak 
AXP J1810-197. 
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During periods of high activity, when bursts occur frequently, the growth of the magnetospheric 
twist may be faster than its decay. The twist may grow to the point of a global instability. When 
the magnetic lines achieve A(/> ^ 1 the magnetosphere becomes unstable and relaxes to a smaller 
twist angle. A huge release of energy must then occur, producing a giant flare; a model of how 
this can happen is dicussed by Lyutikov (2006). It is possible that the 27 December 2004 giant 
flare from SGR 1806-20 was produced in this way — a prolonged period of unusually high activity 
preceded the flare. If the flare indeed released the external twist energy, then the current through 
the magnetosphere should have decreased after the flare, and the coronal luminosity now should 
be smaller than it was prior to the flare. 

A possible way to probe the twist evolution is to measure the history of spindown rate P of 
the magnetar. The existing data confirm the theoretical expectations: the spin-down of SGRs was 
observed to accelerate months to years following periods of activity (Kouveliotou et al. 1998, 1999; 
Woods et al. 2002) . A similar effect was also observed following more gradual, sub-Eddington flux 
enhancement in AXPs (Gavriil & Kaspi 2004). This "hysteresis" behavior of the spindown rate 
may represent the delay with which the twist is spreading to the outer magnetosphere. 

Since the density of the plasma corona is proportional to the current that flows through it, 
changes in the twist also lead to changes in the coronal opacity. The opacity should increase 
following bursts of activity and affect the X-ray pulse profile and spectrum. Such changes, which 
persist for months to years, have been observed following large X-ray outbursts from two SGRs 
(Woods et al. 2001). 

8.7. Further Developments 

1. We focused in this paper on the inner magnetosphere, that is, on field lines that do not 
extend beyond a radius i?max ~ 2i?NS- This part of the magnetosphere is expected to retain most of 
the toroidal field energy, and to dominate the non-thermal output of the star. The outer corona may 
be qualitatively different because it is subject to a strong "cyclotron drag" : the resonantly scattered 
X-rays can escape (rather than convert to e^), thus providing an efficient sink of energy. This drag 
will cause a local change of the electric field in the circuit (TLK; Thompson &: Beloborodov 2005). 
At yet larger radii, beyond ~ 200 km from the star, ions may populate higher Landau levels n > 0. 
Then a new degree of freedom appears in the problem: the transverse motion of the ions. 

2. Like the Solar transition region, the transition layer at the base of a magnetar corona may 
respond to the coronal dissipation by increasing the supply of plasma. In particular, pair creation 
in the transition layer needs further investigation because pairs may feed the coronal current (§ 6.3). 
The state of this layer and its effect on the corona is determined by the level of plasma turbulence 
(see §§ 6.2 and 6.4). This layer may partially insulate the star from the hot corona because plasma 
turbulence suppresses thermal conduction. By trapping the generated heat, the layer may overheat 
and expand to the corona. 
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3. The transition layer also plays a crucial role in the energy balance of the corona because 

it produces radiation that can escape the magnetar. A significant part of the released energy is 
emitted in the hard X-ray band. The rest of the released energy is deposited in an optically thick 
surface layer and reprocessed into blackbody radiaton with ksT < 1 keV. An accurate calculation 
of the ratio of the two emission components (hard X-ray to blackbody) that are fed by the corona 
remains to be done. 

4. The role of transverse waves in the corona, which were suppressed in our numerical exper- 
iment, will need to be studied. The waves produce a three-dimensional pattern of magnetic-field 
fluctuations S{V x B) ^0 and their study will require 3-D or 2-D simulations. 

5. Our plasma model of the corona may be generalized to include the effects of rotation. This 
can be done by changing the neutrality condition p = to p = pco- Note that pco niay be easily 
maintained by a small polarization of the current-carrying plasma because pco ^ c^-c (see eq. [2]). 
It can be taken into account by adding an extra source term 47rpco to the Poisson equation (cf. 
eq. 21), d£^||/d/ = 47r(p — pco)- Then E^y vanishes if p = pco along the magnetic line, so that pco is 
the new equilibrium charge density. 

6. It is important to understand the resistive redistribution of the twist in the magnctosphcre, 
in particular its spreading toward the open field lines, because it is a promising source of torque vari- 
ability in active magnetars (TLK) . This would require the construction of a global, time-dependent 
model of the currents in the magnetosphere. 
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time 




(a) (b) 

Fig. 5. — Spacetime diagrams illustrating the critical character of the breakdown and formation 
of avalanches. The toy model used in this illustration is described in the text; e$e = ^EL = llrrieC^ 
is assumed. Worldlines of particles are shown, starting with one seed electron at the cathode, (a) 
ttres/L = 0.35 (supercritical case), (b) Ures/L = 0.6 (critical case). The magnetar corona is 
expected to evolve close to the threshold for triggering avalanches, when arcs is comparable to 
but somewhat smaller than L (or, equivalently, when the voltage <J>e is just above the minimum 
needed to trigger creation). The avalanches are then initiated stochastically in the regime of 
self-organized criticality. This expectation is confirmed by the numerical experiment described in 
S 5.4. 
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Fig. 6. — Voltage vs. current for circuits with (red circles) and without (black circles) e creation. 
Model A of pair production is represented by filled red circles and Model B by open red circles; 
both models have the resonant Lorentz factor 7res = 10. Each point in the figure represents the 
result of one experiment and shows the time-averaged voltage ^e- The current is expressed in units 
of (see the text). 
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Fig. 7. — Relaxation of the circuit during the first f Hght-crossing times. The circuit parameters 
are: rrii = lOmg, j = 10^ j*, 7rcs = mi/me- Model A is used for pair production (infinite window 
above 7res with free path l^es = 0). A quasi-steady fluctuating state is established after a few L/c. 
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Fig. 8. — Same as Fig. 7, but Model B is used for pair production: narrow window 7res < 7e < 27r, 
with free path ^res = L/2>. 
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Fig. 9. — Circuits with mj/m-e = 10, 30 and various values of 7res ranging from rrii/bme to bmi/rrie. 
Here is the time-averaged voltage and 2iV+ is the time-averaged rate of creation in the circuit. 
Model A of pair production is shown by open symbols, and Model B by filled symbols. 
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Fig. 10. — Fraction of the electric current carried by ions vs. voltage, for the same set of models as 
in Fig. 9. Colors have the same meanings. In the left panel voltage is taken in units of rriiC^, and 
in the right panel in units of 7res"^eC^- 
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Fig. 11. — Snapshot of the circuit with ^^es = rm/me = 10 and Jb = W^j* {Xo/L = 10~^). 
Model A is assumed for creation. The three curves show particle density n = rii + rie + Up 
(dotted), charge density p = e{ni — rie + Up) (solid), and electric field E (long-dash). 
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Fig. 12. — Mean velocities of ions, electrons, and positrons as a function of position in the circuit. 
The mean velocities vx are related to densities nx of the three species X = i, e,p hy vx = jx/vx 
where jx is the current carried by species X. The snapshot is taken of the same model and the 
same time as Fig. 11. 
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Fig. 13. — Dispersion of momentum (around the mean bulk momentum) for ions, electrons, and 
positrons as a function of position in the circuit. The snapshot is taken of the same model at the 
same time as in Fig. 11. Only the corona is shown {h < z < L — h) and the surface layers of height 
h = 0.04L are excluded from the figure. 
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Fig. 14. — Distribution of pair-creation events at 10 < ct/L < 20 over z. The circuit model is the 
same as in Figs. 7 and 11-13. The figure shows separately pair creation by accelerated e~ and 
(solid and long-dashed curves, respectively). 



